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ABSTRACT 


The  objective  of  this  research  was  to  develop  a  new  technique  for  evaluation  of  the 
eigenvalues  of  the  Graetz  problem  in  slip-flow — a  heat  transfer  problem  for  gases  at  low 
pressures  or  in  extremely  small  geometries.  In  this  investigation,  the  velocity  distribution 
with  slip-flow  has  been  obtained,  expressed  simply  in  terms  of  Knudsen  {Kn)  numbers. 
The  expression  shows  that  the  velocity  always  increases  as  the  Knudsen  number  in¬ 
creases.  The  relationship  of  Kn  and  molecular  mean  free  path  for  a  gas  shows  that  Kn  may 
become  large  enough  to  significantly  affect  the  velocity  distribution  and  consequently  af¬ 
fect  the  heat  transfer  properties.  A  mathematical  model  of  temperature  distribution  was 
established  by  combining  the  energy  and  momentum  equations.  A  series  solution  was 
obtained  by  the  method  of  Frobenius.  Also,  expressions  for  the  local  and  overall  Nusselt 
numbers  were  derived.  All  these  expressions  can  be  taken  as  functions  of  Knudsen  num¬ 
bers  and  Graetz  numbers. 

A  new  technique  for  evaluation  of  eigenvalues  for  the  solution  of  the  Graetz  problem 
in  shp-flow  was  developed.  This  method  was  based  on  the  construction  of  a  matrix.  The 
computational  results  show  that  it  is  an  effective  method,  and  the  lowest  five  values  were 
found  for  Kn  from  0.02  to  0. 12.  For  practical  calculations,  relationships  between  eigen- 
values  and  Knudsen  numbers  were  obtained. 
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NOMENCLATURE 


Aw  tube  surface  area  [m^] 

Ok  coefficient  in  Eq.  (3.12) 
bi^k  coefficient  in  Eq.  (4.8) 
c  characteristic  velocity 
Ca  acoustic  velocity 
Cn  coefficient  in  Eq.  ( 1 .3) 
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F  specular  reflection  coefficient 

(Ur  -  Uj)  /  (Uw  -  Ui) 

g  magnification  ratio  in  Eq.  (5.1) 

G  G(r*):  Graetz  function 

Gz  Graetz  number,  {RePr{D/L)) 
hx  local  convective  heat  transfer  coef¬ 
ficient  [W/ra^-K] 

he  average  convective  heat  transfer 
H  characteristic  dimension 
coefficient  [W/m^-K] 


k  thermal  conductivity  [W/m-K]; 

number  of  terms  in  Equation  (5.2) 
Kn  Knudsen  number,  {X/D) 

L  length  of  tube[m] 

Nu  overall  heat  transfer  coefficient, 

(hcD/k) 

Nux  local  heat  transfer  coefficient, 

(hxD/k) 

p  fluid  pressure  [Pa] 
q  heat  flux  per  unit  wall  area  [W/m^] 
Q  heat  transfer  rate  [W] 
r  radius[m] 

r*  dimensionless  radius,  {r/R) 

R  tube  radius[m] 

Re  Reynolds  number,  (guD/p) 

Re*  Reynolds  number  at  Ca 
Rg  ideal  gas  constant  [J/kg-K] 

Pr  Prandtl  number,  (v/a) 

T  T{r,x),  temperature  [K] 
ffi  bulk  temperature  [K] 

7l  temperature  at  x  =  L  [K] 

7)ln,  log-mean-temperature 


difference  (LMTD)  [K] 
u  velocity  in  x  direction  [m/s] 

Mi  average  strearawise  velocity  of  the 
incident  molecules 

Ur  average  stream  wise  velocity  of  the 
reflected  molecules 

f/w  average  streamwise  velocity  of  the 
surface 

V  radial  velocity 

X  distance  along  tube[m] 

X*  dimensionless  distance,  {x/L) 

Greek  Symbols 

a  fluid  thermal  diffusivity,  (k/gc) 
[m^/s] 

Os  thermal  accommodation  coefficient 
(1+4  Kn)  coefficient  in  Eq .  (4. 1 ) 
y  ratio  of  specific  heats 

4l  A  =i-k  in  Eq.(4.8); 

A=(Nu  -  Nux)  in  Fig.  5.9 
r]  mean  free  path  of  gas 


A  eigenvalue; 

X'  eigenvalue  divided  by  g 
fi  dynamic  viscosity[kg/m  s] 

V  kinematic  viscosity  [ra^/s] 

Q  density[kg/m^] 
n  3.141592654 

6  (T-Ty,)fiTo-Tw)  dimensionless  tem¬ 
perature 

iTB-Tw)l(Ta-Tv,)  dimensionless 
bulk  temperature 

6bx  dimensionless  fluid  bulk  temperture 
atx  =  L 
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Subscripts 
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CHAPTER  1 


INTRODUCTION 

1.1  The  Graetz  Problem 

By  the  end  of  the  last  century,  the  problem  of  forced  convection  heat  transfer  in  a 
circular  tube  in  laminar  flow  gained  interest  because  of  its  fundamental  importance  in 
physical  problems  such  as  the  analysis  and  design  of  heat  exchangers. 

The  Graetz  problem  is  a  simplified  case  of  the  problem  of  forced  convection  heat 
transfer  in  a  circular  tube  with  laminar  flow.  With  the  assumptions  of  steady,  incom¬ 
pressible  and  fully-established  flow,  constant  fluid  properties,  no  ’’swirl”  component  of 
velocity,  a  fuDy  developed  temperature  profile,  and  negligible  energy  dissipation  effects, 
Graetz  (1883)  onginally  solved  this  problem  analytically.  The  solution  by  Graetz  in¬ 
volved  an  infinite  number  of  eigenvalues,  and  in  his  paper  only  the  first  two  eigenvalues 
were  evaluated. 

Since  the  accuracy  of  the  Graetz  solution  mainly  depends  on  the  number  of  eigenva¬ 
lues,  it  is  extremely  important  to  obtain  more  eigenvalues,  as  Tribus  and  Klein  (1953) 
pointed  out  For  seventy  years,  the  research  on  this  problem  focused  mainly  on  finding 
more  eigenvalues.  And  Abramowitz  (1953)  employed  a  fairly  rapidly  converging  series 
solution  of  the  Graetz  equation  in  making  the  calculation  and  found  the  lowest  five  values 
with  much  more  accuracy.  SeUars  et  al.(1956)  extended  the  problem  to  include  a  more 
effective  approximation  technique  for  evaluation  of  the  eigenvalues  of  the  problem;  they 

could  get  any  number  of  eigenvalues  as  needed.  This  work  solved  the  Graetz  problem 
completely. 
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1,2  The  Graetz  Problem  in  Slip-Flow 

Applications  of  raicrostructures  such  as  micro  heat  exchangers  have  led  to  increased 
interest  in  convection  heat  transfer  in  micro  geometries.  Some  experimental  work  has 
been  done,  such  as  the  experimental  investigations  in  microtubes  (Choi  et  al..  1991),  in 
microchannels  (Pfahler  et  al.,  1991),  and  in  micro  heat  pipes  (Petersen  et  al.,  1993). 
Therefore,  appropriate  models  are  needed  to  explain  the  significant  departures  in  the 
micro-scale  experimental  results  from  the  thermofluid  correlations  used  for  convention¬ 
al-sized  geometries.  For  example,  the  measured  heat  transfer  coefficients  in  laminar  flow 
in  small  tubes  exhibits  a  Reynolds  number  dependence,  in  contrast  to  the  conventional 
prediction  for  fully  established  laminar  flow,  in  which  the  Nusselt  number  is  constant 
(Choi  et  al.,  1991).  Also,  an  experimental  investigation  of  fluid  flow  in  extremely  small 
channels  showed  that  there  are  deviations  between  the  Navier-Stokes  predictions  and  the 
experimental  observations  (Pfahler  et  al.,  1991). 

Therefore,  some  effects  and  conditions  that  are  normally  neglected  when  consider¬ 
ing  macro-scale  flow  must  be  taken  into  consideration  in  micro-scale  convection.  One 
of  these  conditions  is  slip-flow  (Flik  et  al.,  1992,  Beskok  and  Kamiadakis,  1992).  It  has 
been  found  that  the  analytical  model  combined  with  slip  flow  conditions  can  fit  the  exper¬ 
imental  data  in  microchannels  with  a  uniform  cross-sectional  area  (Arkilic  et  al.,  1994) 
and  with  a  non-uniform  cross-sectional  area  (Liu  et  al.,  1995). 

Slip-flow  occurs  when  gases  are  at  low  pressures  or  for  flow  in  extremely  small  pas¬ 
sages.  At  low  pressures,  with  correspondingly  low  densities,  the  molecular  mean  free 
path  becomes  comparable  with  the  body  dimensions,  and  then  the  effect  of  the  molecular 
structure  becomes  a  factor  in  flow  and  heat  transfer  mechanisms(Eckert  and  Drake,  1 972). 

The  relative  importance  of  effects  due  to  the  rarefaction  of  a  gas  can  be  indicated  by 
the  Knudsen  number,  a  ratio  of  the  magnitude  of  the  mean  free  molecular  path  in  the  gas 
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to  the  characteristic  dimension  in  the  flow  field.  The  effects  of  rarefaction  phenomena  on 
flow  and  heat  transfer  becomes  important  when  the  Knudsen  number  can  no  longer  be 
neglected. 

In  defining  when  slip-flow  occurs,  Beskok  and  Kamiadakis  (1992)  have  proposed 
to  classify  four  flow  regimes  for  gases,  as  follows: 

Continuum  flow:  An  <  lO”^ 

SUp-flow:  10-3  <AGn<  0.1 

Transition  flow:  0. 1  <  ATn  <  1 0 

Free  molecular  flow  10  ^  AT/i 

When  slip-flow  occurs,  the  gas  adjacent  to  the  surface,  in  contrast  to  its  behavior  in 
continuum  flow,  no  longer  reaches  the  velocity  or  temperature  of  the  surface.  The  gas 
at  the  surface  has  a  tangential  velocity,  and  it  slips  along  the  surface.  The  temperature 
of  the  gas  at  the  surface  is  finitely  different  from  the  temperature  of  the  surface,  and  there 
is  a  jump  in  temperature  between  the  surface  and  the  adjacent  gas.  Eckert  and  Drake 
(1972)  give  expressions  for  the  temperature  jump  condition  and  slip  velocity  at  the  sur¬ 
face.  The  slip  velocity  as  a  function  of  the  velocity  gradient  near  the  wall  can  be  expressed 
as  follows: 

^  W  (1.1) 


and  Arkilic  et  al.  (1994)  give  the  expression  as  follows: 


=  IlzE. 
F 


Vr,  (  \ 

^  drjH 


(1.2a) 


(1.2b) 
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which  includes  the  consideration  of  three  accommodation  coefficients  represented  by  the 
specular  reflection  coefficient  F.  For  most  engineering  surfaces,  F  has  values  near  unity. 
In  the  case  of  F  having  a  value  of  one,  Eq.  ( 1 .2)  becomes  Eq.  ( 1 . 1 ).  For  simplicity,  in  this 
investigation,  Eq.  (1.1)  is  applied  to  evaluate  the  velocity. 

The  original  solution  by  Graetz  (which  was  discussed  above)  is  valid  for  continuum 
flow;  however,  for  gases  at  low  pressures  or  in  extremely  small  tubes,  the  flow  may  enter 
the  sli{>-flow  regime,  in  which  case  the  velocity  at  the  tube  surface  is  not  zero.  In  this 
case,  the  heat  transfer  coefficient  depends  not  only  on  the  Reynolds  number  and  Prandtl 
number,  but  also  on  the  Knudsen  number.  This  fact  will  no  doubt  make  the  model  more 
complex  and  the  evaluation  of  its  eigenvalues  more  difficult.  Therefore,  a  new  technique 
is  needed  to  evaluate  the  eigenvalues  for  a  solution  to  the  problem  in  slip-flow. 

1.3  Related  Research 

Graetz  (1883)  originally  solved  the  problem  of  forced  convection  heat  transfer  in  a 
circular  tube  in  laminar  flow,  with  a  developing  temperature  profile.  Figure  1.1  shows 
the  geometry  and  conditions  for  this  problem. 


Fig.  1.1  Coordinate  system  and  thermal  boundary  conditions 
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The  mathematical  statement  of  this  problem  is  as  follows: 

T  =  T{x,r) 


uoC  ^ 


u  =  2M;;,[l-(r/ro)^] 


with  the  boundary  conditions: 

T(R,  X)  = 

T(n  0)  =  To 
T(0,  0)  =  To 

The  nondimensional  form  of  the  problem  is: 

§e_  ^  _j _ 1  a 

dx*  l-r*2r*dr*^  dx*’ 

and  the  boundary  conditions  are  as  follows: 

6(1,  X*)  =  0 
e(r*0)=  1 


for x>  0 
fOTX<0 


for  x<0 


for;c*  >  0 
for  jr*  <  0 


The  solution  for  this  system  can  be  obtained  (Graetz): 


d(r*,x*)  =  I!  C„G„{r*)e-^’‘ 

71  =  1 


where  the  Xi,  are  the  eigenvalues  required  to  make  the  solution  satisfy  the  following  differ¬ 
ential  equation : 

r*Gn"  +  Gn' +X„^r*(\-r*^)Gn  =  0  (1.6) 

Graetz  posed  an  eigenfunction: 


G„(A„)  ^  —  1  -b  -f  + 


(1.7) 
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and  he  gave  only  two  values:  X,i  —  2.704  and  X2  =  6.50.  Unfortunately,  it  was  very  diffi¬ 
cult  at  that  time  to  calculate  the  larger  values  of  Xn 

Sellars  et  al.(1956)  extended  the  problem  to  include  a  more  effective  approximation 
technique  for  evaluation  of  the  eigenvalues  of  the  problem.  They  developed  an  approxi¬ 
mate  method  by  using  three  expressions  to  represent  the  Graetz  functions  in  three  ranges: 
(1)  near  the  center;  (2)  between  the  centerline  and  the  waU;  and  (3)  near  the  wall.  They 
obtained  an  approximate  expression  as  follows: 

Xn  =  4(n-1  )  +  8/3  n=  1,2,3, .  (1.8) 

The  comparison  of  the  values  with  other  investigations  (see  Table  5.3)  shows  that  the 
approximate  method  is  correct  and  effective,  especially  for  larger  n.  The  accuracy,  except 
the  first  eigenvalue,  is  acceptable.  This  work  solved  the  Graetz  problem  completely. 

The  objective  of  this  research  is  to  evaluate  the  eigenvalues  for  the  Graetz  problem 
in  slip-flow  a  heat  transfer  problem  for  gases  at  low  pressures  or  in  extremely  small 
tubes.  To  do  this,  the  velocity  profile  with  slip-flow  must  be  found  first,  and  a  mathemati¬ 
cal  model  of  temperature  distribution  in  slip  flow  must  be  established  by  combining  the 
energy  and  momentum  equations.  Next,  by  using  the  method  of  Frobenius,  a  series  solu¬ 
tion  must  be  obtained  and  finally,  a  technique  for  evaluation  of  the  eigenvalues  for  the 
senes  solution  must  be  developed.  For  practical  calculations,  relationships  between  the 
eigenvalues  and  the  Knudsen  number  should  be  obtained. 


CHAPTER  2 


VELOCITY  AND  TEMPERATURE  DISTRIBUTIONS 

In  order  to  buUd  the  mathematical  model  for  the  Graetz  Problem  in  slip-flow,  the 
velocity  profile  must  be  found  first.  In  this  chapter,  based  on  some  assumptions,  the  ex¬ 
pression  for  velocity  will  be  derived  from  the  continuity  equation  and  momentum  equa¬ 
tion.  The  slip  condition  will  be  used  to  evaluate  the  slip  velocity  and  the  velocity  will  be 
expressed  in  terms  of  Knudsen  numbers.  A  mathematical  model  of  temperature  distribu¬ 
tion  in  slip-flow  will  be  established  by  combining  the  energy  and  momentum  equations. 

2.1  Velocity  Distribution 

As  a  model,  one  can  consider  the  flow  of  a  fluid  in  a  circular  tube  of  radius  R,  shown 
in  Figure  2. 1 : 


L 


Fig.  2.1  Coordinate  system  for  the  problem 


For  this  model  the  following  conditions  have  been  assumed  (Barron,  1994); 

(1)  The  flow  is  steady.  This  means  that  the  properties  of  the  flow  are  time  indepen¬ 
dent 

(2)  The  fluid  is  incompressible  (or,  if  a  gas  is  considered,  the  Mach  number  is  low). 
In  this  case,  the  density  may  be  assumed  constant. 
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(3)  The  flow  is  fully  established.  In  this  case,  the  axial  velocity,  u,  is  a  function  of 
the  radial  coordinate  only,  and  not  a  function  of  the  axial  coordinate.  In  addition, 
the  radial  velocity  is  zero. 

(4)  The  ’’swirl”  component  of  velocity  is  identically  zero.  This  fact  means  that  the 
flow  properties  are  independent  of  the  angular  coordinate  in  cylindrical  coordi¬ 
nates. 

(5)  Fluid  properties  arc  constant. 

(6)  Energy  dissipation  effects  are  negligible. 

(7)  The  tube  wall  temperature  is  constant. 


The  general  continuity  equation  can  be  written  in  cylindrical  coordinates  as  follows: 


§  +  7i:(en')  +  -^(gu)  =  0 


For  steady  flow  of  an  incompressible  fluid,  Eq.  (2.1)  reduces  to: 


For  fully  developed  flow. 


Therefore, 


r  V  =  constant 
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Since  the  radial  velocity  is  zero  at  the  wall  (the  wall  is  impermeable),  we  must  conclude 
that: 

V  =  0  (identically). 

2.1.2  Velocity  Distrihiition  with 
Slip  Condition 

The  Momentum  Equation  can  be  written  in  cylindrical  coordinates  as  follows  (Kays 
et  al.,  1993): 


Qu^  +  Qv^  +  ^  = 

dx  dr  ^  dr 


For  fully-established,  steady-flow  of  an  incompressible  fluid,  the  Navier-Stokes  i 
tion  for  the  axial  direction  reduces  to: 

+  iLA(rd]£)  =  0 


In  this  case,  the  velocity  u  is  independent  of  x,  so  we  can  define  the  constant. 


r  = _ 

‘  AjU  dx 


Then,  Eq.  (2.3)  can  be  written  in  the  following  form: 


)  =  0 


This  expression  can  be  solved  by  directly  integrating  twice  to  yield: 


M  =  C2  /n  r  +  C3  -  Cl 


(2.6) 
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Since  the  velocity  u  is  finite  at  the  center  of  the  tube  (r=0),  we  must  have  C2  =  0. 
At  the  centerline  of  the  tube  ir=0),  the  velocity  is  equal  to  Uc  (the  center-line  veloc¬ 
ity).  Using  this  condition  in  Eq.  (2.6),  we  obtain:  u(0)  =  Uc  =  Cj 

At  the  surface  of  the  tube  {r=R),  the  velocity  is  not  zero  in  slip-flow,  but  is  equal  to 
a  finite  velocity  Using  this  condition,  we  find,  from  Eq.  (2.6): 


U(R)  =  Us  =  Uc-  C]/?2 
Rearranging  and  solving  for  Ci: 

r  =  ~  _ _ 

Making  these  substitutions,  we  find  the  following  expression  for  the 
tribution: 


(2.7) 


(2.8) 

velocity  dis- 


u  =  Uc-(uc-Us)  (r/R)^  =Uc[l-  (r/R)'^  ]  -1-  Us  (r/Rp-  (2.9) 

This  velocity  profile  is  shown  in  Figure  2.2. 

.  r 


Fig.  2.2  Velocity  distribution 


Let  us  now  calculate  the  mean  fluid  velocity  Um  in  the  tube.  The  volumetric  flow  rate 


can  be  written  as  follows: 
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jiF}um  =  I  Qjzrudr 


If  we  introduce  the  dimensionless  variable,  r*,  we  can  write  Eq.  (2.10)  as: 


(2.10) 


Um  =  2j  ur*dr*==2^  [uc -{.Uc-Us)r*\*  dr*  (2.11) 


Carrying  out  the  integration,  we  obtain: 


Mm  =  2[  {Uc/2)  r*^-{ucr-us)  (r*V4)  ] 


Mm  ~  2[  /2  (  Uq  Ug)  /  4  ]  —  f  +  Mj)  /  2 


(2.12) 


The  centerline  velocity  can  be  written  in  terms  of  the  mean  velocity  and  the  slip  velocity. 


as  follows: 


Uc  —  2  U/fi  —  Us 


(2.13) 


Making  this  substitution,  we  find  the  following  expression  for  the  velocity  profile  in  slip- 


flow: 


U  —  2(  Uffi  —  Uy)  (1  —  )  +  Uy 


(2.14) 


If  the  slip  velocity  is  zero,  then  Eq.  (2.14)  reduces  to  the  Poiseuille  distribution; 


«  =  2n^  (1-r^) 


(2.15) 
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lation  of  the  SIid  Velnrif 


The  slip  velocity  can  be  evaluated  from  Eq.  (1.1)  as; 


,,  =  _A/'  \  -  _i  ^MUm-Us)  _ 

=  ^ ^ - 


(2.16) 


Introducing  the  Knudsen  number,  we  find: 


«£  ^  SKn 
««  1  +  SKn 


(2.17) 


Making  this  substitution  in  Eq.  (2. 14),  we  obtain  the  following  expression  for  the  velocity 
profile  in  slip  flow; 


u  _2{l-r*2)  +  SKn 
1  +  SKn 


(2.18) 


The  molecular  mean  free  path  for  a  gas  can  be  calculated  from  the  following  expres- 
sion(Sreekanth,  1968): 


X=^(  -JL-  )i/2 
Q  ^  2RgT  ’ 


(2.19) 


The  acoustic  velocity  for  a  gas  is  given  by: 


=  (y/?57)i/2 

where:  y  =  c/cy.  Substituting  Eq.  (2.20)  into  Eq.(2.19),  we  obtain: 


(2.20) 


2  =  JLr  ^  ^l/2  _  1.4829/t 
QCa^  2  ’  QCa 


(2.21) 
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for  a  gas  having  a  specific  heatratio,  y  —  1.40.  The  Knudsen  number  may  then  be  written 
as  the  reciprocal  of  the  Reynolds  number  based  on  the  sonic  velocity  in  the  gas  as  follows: 


Kn  = 


^  ^1/2  =  _I_ 
Dgca^  2  ^  Re* 


( f  = 


0.6743  Re  * 


(2.22) 


From  the  expression,  it  is  obvious  that  the  Knudsen  number  is  dependent  on  Re  *  for  a  gas. 
Thus,  for  a  given  temperature,  Kn  can  be  calculated  from  Eq.  (2.22)  and  the  velocity  dis¬ 
tribution  can  be  determined  from  Eq.  (2.18). 

For  example,  for  nitrogen  gas  at  300  K  (26.8°C  or  80°F)  and  atmospheric  pressure, 
the  gas  mean  free  path  may  be  determined  from  Eq.  (2.21),  where  the  property  values  for 
nitrogen  gas  are  as  follows: 

0  =  1.6332  kg/m3  =  0.07106  Ibm/ft^ 
p  =0.01784  mPa-s  =  0.04316  Ibm  /ft-hr 
Rg  =  206.8  J/kg-K  =  55.15  ft-lbfrlbm-°R 


A  =  M1784)  (10-3)  _  ^ 

(1.6332)  ^  (2)  (206.8)  (300)  ^ 


A  =  0.2189  X  10“6  m  =  0.2189  pm 


If  we  accept  a  difference  of  5  percent  between  the  case  for  slip  flow  and  the  case  for 
continuum  flow  as  the  effect  of  the  slip  condition,  from  Eq.  (2. 17)  we  find  that  the  slip 
flow  effects  become  significant  when: 


«£  ^  SKn 
1  -f  SKn 


0.05 


Then 


14 


%Kn  ^  0.05  ,  or  Kn^  0.00625 
The  corresponding  tube  diameter  is: 

D  =  0.2189/0.00625  =  35.0  pro 

For  slip-flow  (10~^  <  Kn  <  0.1),  the  tube  diameter  range  is: 

^min  ~  2.2  pm  for  Kfi—0. 1 

^max  =  218.9  pm  forA>i=0.001 

and  from  Eq.  (2.18)  the  maximum  velocity  can  be  found  in  the  range  of: 
l^niax  /  Um  =  1.556  for  Kn  =  0.1 

t^max !  Uju  =  1.992  for  Kn  —  0.001 
while  Umax  /  Urn  =  2  for  Kn  =  0  (no  slip) 

which  means  that  the  velocity  difference  between  the  wall  and  the  centerline  can  be  re¬ 
duced  significantly  in  small  size  tubes. 

2.2  Temperature  Distribution 

The  general  temperature  field  equation  for  flow  of  an  incompressible  fluid  with  zero 
swirl  or  angular  components,  zero  energy  generation,  and  negUgible  frictional  energy  dis¬ 
sipation  is  as  follows: 
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dt  dr 


wil)  =  +  £li 


(2.23) 


For  steady  flow  and  for  zero  radial  velocity  (v  =  0),  the  energy  equation  reduces 
to: 


u—  =  a  r  - 
dr  “  ^  r 


A( 

dr^  dr 


-gi 


(2.24) 


From  an  order-of-raagnitude  analysis  for  the  case  in  which  the  tube  length  is  much  larger 
than  the  tube  diameter,  the  second  term  on  the  right  side  of  Eq.  (2.24)  is  much  smaUer 

than  the  first  term  on  the  right  side  and  the  energy  equation  can  then  be  written  in  the  fol¬ 
lowing  form: 


r  dr 


(2.25) 


One  can  then  consider  the  case  for  flow  which  is  fully  established  hydrodynamically 
at  the  end  of  the  insulated  section  ( a  long  entry  length),  but  which  is  developing  thermally 
due  to  the  temperature  jump  on  the  tube  wall.  The  physical  model  is  illustrated  in  Figure 
1.1.  The  velocity  profile  is  fully-established  at  the  end  of  the  insulated  section,  x=0,  and 
the  temperature  of  the  fluid  entering  the  uninsulated  section  is  uniform  T  =  %  .  The 
boundary  conditions  for  this  situation  are: 


T(R,x)  =  T^v  forx>0 


T(r,0)=To  forx^O 

Using  the  dimensionless  variables,  the  energy  equation  may  be  written  in  the  following 
form: 
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dO  _  1  a  ,  *86  X 

aL  dx*  r  *  dr  *  .  dr*  ’ 


(2.26) 


Using  the  velocity  distribution  given  by  Eq.  (2.14) 

U  =  2(  Unt  -  Us)  (]  -r*^ )  +  Us  =  2(  Um  -  Us)  (1  -r*  ^  +  4Kn)  (2.27) 

the  energy  equation  reduces  to: 

UrrX^^  BB  _  Id  B6 

2aL{  1  +  %Kn)  dx*  ~  (l-r*^  +  4Kn)r*dr*^  ^*Tr*  ^  (2.28) 

The  Graetz  number  Gz  is  defined  by: 

Gz  =  Re  Pr  (D/L)  = 

P  k  L  aL 

The  energy  equation,  Eq.  (2.28),  can  be  written  as  follows: 

I  Gz  dd  _  _ _ 1  d  /  dd 

2(  1  +  8^71)  dr*  (l-r*2  +  4A:n)r*dr*  ^  ^  dr*  ^  (2.29) 

with  boundary  conditions: 

6(l,x*)  =  0  forx*>0 

B(r*0)=:]  forx*<0 


2.3  Summary 

In  this  chapter,  the  velocity  distribution  with  slip-flow  has  been  obtained.  It  can  be 
expressed  simply  in  terms  of  the  Knudsen  number.  From  the  expression,  it  is  obvious  that 
the  velocity  increases  always  as  Kn  increases.  From  the  relationship  of  Kn  in  term  of  mo¬ 
lecular  mean  free  path  for  a  gas  A  and  diameter  of  tube  D,  we  can  see  that  Kn  in  microtubes 
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may  become  large  enough  to  significantly  affect  the  velocity  distribution  and  consequent¬ 
ly  affect  the  heat  transfer  for  this  problem.  Also,  a  mathematical  model  of  temperature 
distribution  in  slip-flow  has  been  established  by  combining  the  energy  and  momentum 
equations. 


CHAPTER  3 


ANALYTICAL  SOLUTION 

In  the  last  chapter,  the  velocity  distribution  was  expressed  in  terms  of  mean  velocity 
and  Knudsen  number,  and  a  mathematical  model  of  temperature  distribution  in  slip-flow 
was  established  by  combining  the  energy  and  momentum  equations.  In  this  chapter,  a 
series  solution  will  be  obtained  by  the  method  of  Frobenius.  Considering  the  given 
boundary  condition,  a  temperature  distribution  in  terms  of  a  generalized  Fourier  series 
will  be  derived.  Also,  expressions  for  the  local  and  overall  Nusselt  numbers  will  be 
obtained. 
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and 
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)+^'<  I-'-*"  +4Kr,)G=  0 


with  boundary  conditions: 

G(l)  =  0, 


G(0)  =  1. 

The  solution  of  Eq.  (3.2)  is; 

The  constant  C  in  Eq.  (3.4)  will  be  evaluated  below,  and  X  will  be  evaluated  in  the  next 
chapter. 

The  solution  of  Eq.  (3.3)  may  be  obtained  by  the  method  of  Frobenius.  Suppose  we 
take  the  function  G(  r*)  as  a  power  series. 


G  (r*)  =  2  Oj  r*^ 


Then, 


y=0^ 


G'  (r*)  =  ^jo:  r*j~^=  ^ 
y=i  ^  j=c 


2,  (1  +  l)a.^j  r*j 


G"  (r*)  —  ^(j  +  l)j  fl.  .  r*j  2 


^(J  +!)(/■  +  2)a.+2  r*j  (3.7) 


Making  these  substitutions  into  Eq.  (3.3),  we  obtain: 

^  (j  +  1)  (/  +  2)aj^2  2  (j  +  1)  fl.  r*j~^  + 

J-'J  ;=0  •' 


+  X^[  (1  +  4Kn  )  S a:  r*j-  2  a-  1  = 

;=o  ^  ■'  ^ 
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Expanding  Eq.  (3.8)  and  multiplying  each  term  by  r*,  we  obtain: 

(2fl2^*  +  Sa^r*^  +  12a4r*^  +  ...)  +  {a^  4-  2a2^*  +  +  Aa^r*^  +  ...)  + 

+  X^[  (1  +  AKn  )  {a^r*  +  +  02^*^  +  ...)  -  +  ...)]  =  0  (3.9) 

The  first  two  constants  oq  and  ai  are  arbitrary,  so  let  us  take  the  even  solution  with 
the  following  values: 

flo  =  1  and  fl]  =  0 

Equating  the  coefficients  on  like  powers  of  r*  in  Eq.  (3.9),  we  obtain: 

2  a2  +  2  a2  +  X^{  1  +  AKn  )  aq  =  0 
or, 

a2  =  -(A/2)2(l  +AKn) 
and 

6  as  +  3  013  +  (  1  +  AKn  )  aj  =  0 

or, 

as  =  0 

In  fact,  we  find  that  all  terras  involving  odd  numbered  subscripts  drop  out. 
^3i2k-i=0  fork  =  1,2,3, ... 


For  the  coefficients  aj,  withy  >4,  we  find  the  following  recursion  relationship; 
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(3.13) 
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Making  these  substitutions  into  Eq.  (3.13),  where  dG(l )/dr*  =  dG(r*)/dr*\r*=  1 ,  we 

obtain  the  following  condition: 

00 

1  +  Say  [  1  +  2;  (  )  ^  ]  =  0  (3.14) 

Eq.  (3.14)  defines  the  present  problem.  The  coefficients  a2j  are  functions  of  the 
eigenvalues  where  n  =1,2,3, ... .  The  eigenfunctions  for  this  problem  can  be  written 
as  follows: 

00 

G„(r*)  =  Say  (A„)  r*y  n  =  1,2,3,...  (3.15) 


3.1.3  Determination  of  CnRstantx 

We  can  write  the  solution  for  the  temperature  distribution  in  terms  of  a  generalized 
Fourier  series,  as  follows: 


e  (r*,x*)  =  S  CnGnir  *)  exp  [  -  —  ^  x*  (1  +  S  Kn) 
n=l  Gz  ^ 


(3.16) 


Note  that  the  lower  limit  on  n  is  now  1,  which  is  arbitrary. 

The  constants  can  be  found  from  the  entrance  condition. 


atA:=0(orj:*=0);  T(  r,  0  )=  Tq,  or  e(r*0)=  1 


Making  this  substitution  into  Eq.  (3.16),  we  obtain: 

00 

SCn  Gn(r*)  =  1  p 

The  governing  differential  equation,  along  with  the  boundary  conditions,  is  a  Sturm- 
Liouville  problem,  with  a  weight  function, 


23 


w  =  r*(  1  -  r*-  +  4  ) 

and  from  orthogonality, 


I  0°" 


Gmf  *  (  i-r  +  ^Kn  )dr  *  =  0  for 


The  constants  may  be  evaluated  by  multiplying  Eq.  (3.17)  on  both  sides  by 


Gmi"  *  (  1  -  r  +  AKn  ) 


and  integrating  between  0  and  1.  Only  the  term  in  which  m  =  n  is  non-zero,  and  we  find: 


G„  r*(  l-r*2  AKn)  dr*=  (Gn)^r*(  1-r*^  +  AKn  )  dr*  (3.18) 


The  integral  on  the  leftside  may  be  evaluated  from  the  differential  equation  (Eq.  (3.3)), 


as  follows: 


Gn  r*(  l-r*2  +  AKn)  = 


_] _ ) 

(A„)2dr*  ''  dr*  ’ 


Gn  r*{  l-r*2  +  AKn  )  dr*  = 


— L.(  )  I  1 

ar*  '  'o 


which  becomes 


Gn  r*i  l-r*^  +  AKn  )  dr*=  - 


1  /  dGn  . 

(A„)2^  dr*  V  =  1 


(3.19) 
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The  integral  on  the  right  side  of  Eq.  (3.18)  can  be  evaluated  from  (Graetz,  1885,  see 
Appendix  D) 


l-r*2  +  4Kn  )  dr*=  :^[  (  |^  )  |^  ] 

2  Afi  Sa  /I  dr 


Making  the  substitutions  for  the  integrals,  we  find  the  following  values  for  the  constants 
in  the  series  expansion: 


Cn  = 


OA  „  I 


(3.20) 


The  term  in  the  denominator  may  be  evaluated,  as  follows: 


[(f)]  r^V)] 

dA  n  r*=  1  J=0  dX  n  1  y=0^  dX  n 


Each  of  the  terms  in  this  equation  may  be  worked  out  from  the  previous  results,  Eq. 
(3.10),  as  follows: 


^2,  _  .A 


(_)^  =  _(_^)  (  1  +  4  ) 


(_rr±.)  = 
^  dX  \ 


(■^)  [1  +  (  1  +  4  )2] 


(-ZIL) 

^  dX  \ 


-  An  (  1  +  A  Kn  ) 


[  5  +  (1  4-  4  ATn  )2] 
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For  the  special  case  of  x/L/Gz  >  0.05,  which  means  that  the  entrance  effect  can  be 
neglected  for  the  rest  of  the  tube,  only  the  first  term  in  Eq.  (3.15)  is  needed  to  represent 
the  temperature  distribution  accurately.  In  this  case,  we  have: 

a  *  „  *\  _  -^)  ~  r  2(Aj)a:*(1  +  8  Kri) 

'  ’  ' - fpr^ - Cfi,(r*)  exp[ - ! - g- - (3.21) 

The  expressions  for  the  temperature  distribution  are  summed  as  follows: 


e  (  r*,x*)  = 


^  CnGnir*)  exp  [ 

n-1 


2  {Xn)^  X*  (1  +  8  Kn) 


(3.16) 


where 


Gn{r*)  =  ^a2j  (Xn)  r*^j 
j=0  ^ 


n  =  1,2,3,  ... 


(3.15) 


C„  =  - 


-*“1  <  It  >  1 

n  T*=  1 


(3.20) 


[(f)]  -2(^) 

OA  n  r*=  \  y=0  UA  n 


From  these  expressions,  we  can  see  that  the  coefficients  ay  and  must  be  predetermined 

in  order  to  calculate  the  temperature. 


26 


3.2  Heat  TVansfer  Coefficient  Correlation 
3.2.1  Bulk  Temperature 

The  bulk  or  average  temperature  can  be  determined  from: 


(«/«„)  6  r  dr 


where: 


=  2- j"  {u/u„ 


i)  0  r*  dr* 


(3.22) 


_ 

-  r — t” 


Using  the  velocity  distribution  from  Eq.  (2. 1 8)  and  the  temperature  distribution  from 
Eq.  (3.16),  we  can  evaluate  the  bulk  temperature  at  any  location  along  the  length  of  the 


tube,  as  follows: 


Bjj  =  ^ -  pyp  r  ^  X*  (1  +8  Kn) 

^  «=i(l  +  8^:n)  ^  Gi - 


Gn  (  1-r  +  dKn  )  r*dr*  (3.23) 


The  integral  had  been  worked  out  previously  in  Eq.  (3.19).  Making  this  substitution, 
obtain  the  expression  for  the  local  bulk  temperature  of  the  fluid  in  the  tube. 


6.= _ 4 _ y  _CiL.  ,  2  (A„)2  X*  (1  +  8  Kn) 

^  (1  +  8/:n)„ti  a„)2  [ - Gi - ^ 


]  (3.24) 
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.2.2  Local  Heat  Transfer  Coefficient 


The  local  or  ’’point”  convective  heat  transfer  coefficient  can  be  defined  by: 
Q/Aw  =  hx  (  Tb  -  Tw  )  (3.25) 


The  heat  flux  can  also  be  written,  as  follows: 


=  — it  (  —  1  I 
Aw  *  ^  dr  ''  =  R 


k  (  Tq-  Tyy  )  J0 
R  ^  dr* 


(  )  I 

^  a;.  *  ^  V  =  1 


(3.26) 


Equating  the  heat  flux  from  Eqs.  (3.25)  and  (3.26),  we  obtain  the  expression  for  the  local 


convective  heat  transfer  coefficient. 


= _ L.  (  )  I  -  2  k  f  de  .  , 

^  Rds^  dr*  =  -Di:  ^  ^  1 


(3.27) 


Making  the  substitutions  from  Eq.  (3.16)  for  the  temperature  gradient  at  the  wall  and  Eq. 
(3.24)  for  the  local  bulk  temperature,  the  following  expression  is  obtained  for  the  local 


or  ’’point”  Nusselt  number. 


Nu,  =  ^  = 


n  =  l  (J7  J 


(3.28) 


(1  +  SKn)  2  C,^%^exp  [  -  ^  0 +8  ATn) 

_  _ _ /I  =  1  ar^  Qz  J 

2  a.r  xy  ±sm 

n  =  i(4J"  dr*  ^  Gz  J 


(3.29) 
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For  the  special  case  of  x*/Gz  >  0.05,  Eq.  (3.29)  reduces  to  the  following 

to 

Nux  =  (  1  /  2  )  (1  +  8  /(/I )  (  )2 


expression: 

(3.30) 


3«2,3  Overall  Convective  Heat  Transfer 
Coefficient 

The  average  or  overall  convective  heat  transfer  coefficient  is  defined  through  the 
following  expression; 


Q  =  h^{jt  D  L)  {A  7  ),M  = 


hx  (  To-r^)  {n  D  dx) 


(3.31) 


where  ( AT )u^  -  log-mean-teraperature  difference  (LMTD). 

The  LMTD  may  be  written  in  terms  of  the  inlet  temperature  Tq  and  the  exit  bulk 
temperature  Tl  ,  as  follows: 


{A  T  )rj.  - - LTo  T'yv  )  -(  Tl-  ^  )  (  Tq-  ) 


/«  [  (  7-0-  )  /  (  7^-  7,  )  ]  - 

Let  us  define  the  dimensionless  LMTD,  as  follows: 


In  (  0g^) 


(3.32) 


G  J  M  “ 


T)Uy 

To-  Tyy. 


(3.33) 


^LN  - 


(  ^B.L~  ^  ) 

) 


(3.34) 


The  expression  for  the  average  convective  heat  transfer  coefficient  can  then  be  written, 
from  Eq.  (3.31): 


hr  9 r,  dx* 


(3.35) 
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We  had  obtained  the  following  result  previously  in  Eq.  (3.28): 


u  _  2  k  f  d6  \  I 

^  ^  r*  =  1 


h, e,  =  -1,^  ^exp  ( -  ^*(1  +  8  m 


(3.36) 


Making  this  substitution  into  Eq.  (3.35)  and  integrating,  we  obtain  the  following  result 
for  the  average  Nusselt  number: 


Nu  =~^ 
k 


Nu  = - ^ _ 

^z^l  +  SKn) 


1:^ 

n=i(A„)2  dr*  ^  Gz  ^ 


(3.37) 


The  expression  for  the  average  Nusselt  number  may  be  written  in  a  somewhat  more 
compact  form,  as  follows.  At  the  inlet  of  the  tube,  the  temperature  of  the  fluid  is  uniform, 

T  (  r,  0 )  -  Tq  ,  6b,  0  =  1. 

Using  Eq.  (3.24),  we  obtain: 


=  1  =  _ 4 _  ^_Cn.  dGn{\) 

^•0  (1  +8/^/z)„fi(A„)2  dr* 


(3.38) 


and. 
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fl  =  4  y  C„  JG„(1) 

(1 +  8i&i)„ri(A„)2  "57^ 


By  comparison  with  the  expression  given  in  Eq.  (3.37)  and  using  Eq.  (3.34),  we  obtain: 


—  Gz  {  \  —  6dj  )  . 

= - 4  0,  '  =  -|  h  (  V  ) 


(3.40) 


The  expressions  for  both  local  and  average  Nusselt  number  are  summed  as  follows: 


{\+^Kn)  I  [  -  AMilli+lM  i 

_  _ _ n  =  l _ Cir  _  J 

2^  ^  [  _  _2a.)^A:»(l+8^ 

n=\{Ky  dr*  Gz  ■' 


(3.29) 


Nu  =  — . .  - y.  SiL.  2  (A„)2  X*  {1+8  Kn\. 

0iA\  +SKn)^i{X„)2  dr*  - Gi - 


(3.37) 


(  ^B,L~  ^  ) 


(3.34) 


4  y  Cn  dGrll) 

(1  +  8/:n)„t,  0^2  ^2cp[ 


2  (A„)2  X*  (1  +  8  Kn) 


(3.39) 


Gn{f  *)  —  ^  ^2/  (^rt)  ^ 


“2/ 

y=0 


«  =  1,2,3, 


(3.15) 
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C„  =  - 


2 


[(f)] 

OA  n  1 


= 

y  =  0^  dX 


) 


(3.20) 


_  dG„{r  *)  _  yn- 

dr*  '^  =  1  “  ^2y 

From  these  expressions  we  can  see  that  the  coefficient  a2;  and  must  be  predetermined 
in  order  to  calculate  the  Nusselt  numbers. 

3.3  Summary 

In  this  chapter,  a  series  solution  for  the  mathematical  model  of  temperature 
distribution  in  slip  flow  in  circular  geometries  has  been  obtained  by  the  method  of 
Frobenius.  Considering  the  given  boundary  condition,  a  temperature  distribution  in 
terms  of  a  generalized  Fourier  series  has  been  derived.  Also,  expressions  for  the  local 
and  overall  Nusselt  numbers  have  been  obtained.  All  these  expressions  can  be  taken  as 
functions  of  the  Knudsen  number  and  the  Graetz  number.  In  order  to  calculate  either  the 
temperature  or  the  Nusselt  numbers,  the  coefficient  azj  and  must  be  predetermined. 


CHAPTER  4 


EVALUATION  OF  EIGENVALUES 

In  the  last  chapter,  we  obtained  a  series  solution  for  the  temperature  distribution. 
Also,  expressions  for  the  local  and  overall  Nusselt  numbers  have  been  obtained,  as  func¬ 
tions  of  the  Knudsen  number  and  the  Graetz  number.  In  this  chapter  we  presents  a  tech¬ 
nique  for  expansion  of  the  coefficient  a2j  and  evaluation  of  eigenvalues  for  the  solution  of 
the  Graetz  Problem  in  slip-flow,  since  the  coefficient  a2j  and  eigenvalues  roust  be  pre¬ 
determined  for  the  calculation  of  either  the  temperature  distribution  or  the  Nusselt  num¬ 
bers.  A  matrix  will  be  constructed  and  a  formulation  described  to  find  the  coefficients  a2j 
du-ectly  as  well  as  4.  Based  on  these  4,  the  eigenvalues  Xn  can  be  calculated  easily. 

4.1.  Introduction 

The  senes  solution  of  Eq.  (3.8)  can  be  expressed  as  Eq.  (3.16),  which  required  the 
solution  of  Eq.  (3.10). 

Letting  p  =  (  \+AKn  ),  Eq.  (3.10  )  can  be  rewritten  as  follows: 

r*^)G  =  0  (4.1) 


with  boundary  conditions: 

G(0)  =  1 
G(l)  =  0 

A  particular  solution  of  Equation  (4.1)  satisfying  condition  (4.2)  is 

00 

G(r  *)  =  2  ajT 
k=o 


(4.2) 

(4.3) 

(4.4) 
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Substituting;  =  2k,  the  recursion  relationship  Eq.(3.10)  can  be  rewritten  as  follows: 


ao  =  \  ,  4ai +PX^=0 

}  (4-5) 

4k\  +  -  ak-2]  =  0 

f0Tk=  2,3.4,... . 

From  Eq.  (4.5)  we  can  see  that  the  expansion  of  ok  involves  A  with  different  powers. 
In  order  to  satisfy  Eq.  (4.3),  the  parameter  A  must  take  on  an  infinite  number  of  discrete 

values  A., .  and  Gf  U  can  be  taken  as  a  fnnction  oM.  Then,  an  eigenfuncdon  is  obtained  by 
rearranging  terms  according  to  A^^- 


G(f)  I...,  -  OO.X)  =  =  1  +  +  d^i’'  +  ...  =0  (4.6) 

One  of  the  procedures  employed  in  determining  the  values  ofX  is  to  (1 )  expand  each 
coefficient  at  to  a  large  k  by  relationship  Eq.  (4.5);  (2)  add  up  all  the  terms  with  the  same 
order  ofA  in  determining  the  coefficients  dt  in  Eq.  (4.6);  and,  (3)  determine  the  valnes  of 
A„.  This  technique  is  computationally  inefficient  for  large  k  because  the  expansion  of  nr 
for  large  k  is  very  complicated  and  makes  the  determination  oU„  from  the  eigenfunction 
C(  1,  A)  =  0  extremely  difficult  It  is  the  purpose  of  this  chapter  to  obtain  a  formulation  for 

calculation  of  nr  in  Eq.  (4.5)  as  well  as  forA„  in  the  eigenfunction  expansion  (Eq.  (4.6))  in 
a  more  effective  way. 

4.2.  Formulation  of  <4 

4^2.1  Expansion  of  a..  fATfi  =  0  nr  p  =  1) 

Firstofall,letusexpandthecoefficientsqr.  Forsiraplification.letcoefficientsni  be 

expanded  from  the  recursion  relationship  Eq.  (4.5)  for/5  =  1  (or  Xn  =  0)  as  follows: 
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<30  =  1 


aj  =  -  X2/[2(2)(1)(1!)2] 


a2=  [X222+X4]/[2(2)(2)(2!)2] 


03  =  -[  X‘^(22+42)+X6]/[2(2)(3)(3!)2 

04  =  [  X42262+X6(22+42+62)+X8]/[2(2)(4)(4!)2j 

as  =  -  {  X^[2M±S!aifc£}l+X8(22+42+62+82)+Xi0}/[2(2)(5)(5!)2] 

06  =  {X62262102+  X8[2262+82(22+42)+io2(22+42+62)]  +X10(22+42+...+io2)+ 

+X12}/[2(2)(6)(5!)2] 

07  =  -{  X8[2262io2+I22[2262  +  82(22+42)]  +X10[2262+82(22+42)+i 02(22+42+62)+ 
+122(22+42+62+82)]  +X12(22+42+.  .+122)+X14}/[2(2)(7)(7!)2j 

08=  {  Xg2262l02l42+XlOr2262l02+122r2262  +  S2r22+42^1+142r2262+ 

±iQ^I2^±4!±fiiai+ 2[2262+82(22+42)+ 1 02(22+42+62)+ 1 22(22+42+62+^2)+ 
+  1 42(22+42+...+ 102)]  +X^4(22+42+__^142)+)^]6}/[2(2)(8)(8|)2] 

.  (4.7) 

The  underlines  are  explained  in  examples  1  and  2  below. 

From  these  expansions  of  we  can  see  the  complexity  of  the  expressions  with  large 
k.  Both  the  number  of  terms  and  the  constants  increase  with  increasing  k.  Because  the 
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denominator  increases  very  rapidly,  it  will  be  extremely  difficult  to  find  the  a^'s  we  need 
for  evaluation  of  eigenvalues  X„ . 

To  see  more  clearly  the  relationship  between  Ok  and  4,  let  us  construct  a  matrix  from 
the  expansions  of  Ok. 

4.2.2  Construction  of  the/>.  r  Mafriv 

Since  we  need  to  add  up  the  numbers  related  to  the  order  of  A^in  ok  the  expansions  of 
ak  can  be  rearranged  as 


k 

1 

2 

3 

4 

5 

6 

7 

i 

a^) 

a^) 

1 

1 

4 

0 

0 

0 

0 

0 

0 

2 

4 

64 

JL 

64 

0 

0 

0 

0 

0 

3 

0 

-20 

2304 

-1 

2304 

0 

0 

0 

0 

4 

0 

144 

56 

I 

147456 

147456 

147456 

0 

0 

0 

5 

0 

0 

-1424 

-120 

-I 

14745600 

14745600 

14745600 

0 

0 

6 

0 

0 

14400 

7024 

220 

1 

2123366400 

2123366400 

2123366400 

2123366400 

0 

7 

0 

0 

0 

-219456 

-24304 

-364 

-] 

416579814400 

416579814400 

416579814400 

416579814400 

8 

0 

0 

0 

2822400 

1596160 

67424 

560 

106542032486400 

106542032486400  106542032486400 

106542032486400 

so  a  matrix  B  is  constructed; 


■^.1 

0 

0 

0 

0 

0  ...' 

^2,1 

^2,2 

0 

0 

0 

0  ... 

0 

^3,2 

^3.3 

0 

0 

0  ... 

0 

^4,2 

^4,3 

^4,4 

0 

0  ... 

0 

0 

^5,3 

^5.4 

^5,5 

0  ... 

0 

0 

^6,3 

^6.4 

^6.5 

^6.6  - 

... 

•  «  * 

•  •• 

**• 

From  the  matrix  and  the  expansions  a^,  we  observe  that 


b  = 

1.1  2(2)(1)(1!)2 

b  = 

2.2  2(2)(2)(2!)2 


^2.1  -  ^2,2  2^  =  ^22  =  Z?22  (5  )2 

^,  =  1 

b  = 

3.3  212)(3)(3!)2 

h,2  =  ^3.3  (2"  +  42)  =  ^>33  2(2X3-2)(i2  +  22)  =  2^2X1)  j;  ^2 

5,  =  1 

b  = 

4.4  2(2)(4)(4I)2 

^4.3  =  ^4.4  (2^'  +  42  +  62)  =  b^^  2^2)(4-3)^22  ^  22  +  32-) 

=  ^,4  2(2)(l)(l2  4-  22  +  32)  =  2(2X1)^  (5,)2 

f,  =  l 

*4,2  =  *4.4  <2%')  =  *4,4  20*«)(3^P)  =  2»M>  [i  +  1)2[2'‘  (.  )2]] 

S.^l  J,=l 
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where  we  define  the  summation  for  the  upper  limit  of  zero  as  follows: 


^  (^i)^  =  0 

j,  =  i 

From  these  observations,  we  can  deduce  a  formulation  to  calculate  bi^i,  directly  re¬ 
lated  to  the  index  instead  of  from  the  expansion  of  each  term  The  formulation  is: 


22‘‘(/!)2 


and  for  i>  k 


bi,k  ^  X  (^2  +  1)^[  2" 


(4.8) 

where  A  —  k,  and  then  the  coefficients  in  the  eigenfunction  can  be  found  from 


2k 

‘>k  = 


Example  1: 

^5.3  =  ^5,5  22(5-3)[i;  (^2  +  5j2]  ] 

^2  =  2  ^,  =  1 

~  22(5)(5!)2  ^ 

^  14, 74^  600  +  4^(1^  +  2^)] 


(4.9) 


-1,424 
14, 745,600 
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Example  2: 


(-!)« 

22(8)(g!)2 


(52  +  1)^[  Z  51^]  +6^Z  (52  +  1)2[  Z  5,2]  + 

^2-2  J,  =  l  5^  =  2  j,  =  I 


5-1  ^2-1 

+  1^Z  (52+  1)2[2'5i2]  ] 

•^2  ”  2  1  “  1 

=  Jio^  2^5^32^  +  +  42(12  +  22)]  + 

+  +  42(l2  +  22)  +  52(12  +  22  +  32)]} 

=  1,596,160 

106, 542,032,486,400 


These  two  examples  are  also  shown  as  the  underlined  parts  in  expansion,  Eq.  (4.7); 
when  applying  the  recursion  relationship,  Eq.  (4.5)  and  the  procedures  described  in  part 
4.1,  the  computations  are  not  easy  but  rather  are  complex  and  time-consuming.  For  ex¬ 
ample,  the  computation  of  b8,5  involves  the  expansions  of  a^,  02, ...,  as  and  addition  of  all 
term  with  A  ^0  in  those  expansions.  Thus,  these  examples  also  show  that  the  computation 

by  the  formulation  in  Eq.  (4.8)  is  much  more  effective  than  that  found  by  the  recursion 
relationship,  Eq.  (4.5). 

4.2,3  Expansion  of  gu  fATw  ^ 

From  the  above  two  examples,  we  can  see  the  correctness  and  effectiveness  of  Eq. 
(4.8)  and,  based  on  it,  the  formulation  for  A>i  >  0  also  can  be  derived. 

For  P  >1  the  corresponding  expansion  of  are  as  follows: 
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ao  =  1 


ai  =  -  X2p/[2(2)(l)(i!)2] 

02=  [X222+X4p2j/f2(2)(2)(2!)2] 

03  =  -[  X4p(22+42)+X6p3]/[2(2)(3)(3!)2] 

04  =  [  X42262+X^p2(22+42+52)^.^8p4j/|'2(2)(4)(4 |)2j 

05  =  -  {  >.^P[2262+82(22+42)]+X8p3(22442+62+g2)^.Xl0p5}/f2(2)(5)(5,)2] 

06  =  {  X62262l02+X8p2[2262+82(22+42)+io2(22+42+62)]+xlOp4(22+42+  +102) 

+Xl2p6}/[2(2)(6)(5j)2j 

07  =  -{  X8p[2262l02+122[2262  +  82(22+42)]  +  XJ0p3|-2252+82(22+42)+ 

+102(22+42+62)+122(22+42+62+82)]  +Xl2p5(22442+  +122)+ 

+X>4p7}/p(2)(7)(7j)2] 

08=  {  X82262l02l42+XlOp2p262io2+i22[2262  +  82(22+42)]+ 

+142[2262+  82(22+42)+!  02(22+42+62)]]+  Xl2p4[2262+82(22+42)+ 

+ 1 02(22+42+62)+ 1 22(22+42+62+82)+ 1 42(22+42+^+ 1 02)]+ 
+Xl4p6(22+42+__+i42)+)^16p8}/[2(2)(8)(8!)2] 


(4.10) 


so  the  matrix  is  modified  as  follows: 


B'  = 


0 

0 

0 

0 

0  ...■ 

K\ 

0 

0 

0 

0  ... 

0 

^3.2^ 

0 

0 

0  ... 

0 

^4,2 

^4,3^^ 

0 

0  ... 

0 

0 

^5,3^ 

b5,5^' 

0  ... 

0 

0 

^6,3 

^6,4/32 

- 

•  •  #  •  *  • 

(4.11) 


2k 

*  i=k 


(4.12) 


All  the  expressions  for  evaluation  of  eigenvalues  of  the  Graetz  Problem  in  slip-flow 
are  summed  as  follows: 


G{r*)  1^^,  _  G(1,A)  =  Z  =  1  +  +  d^X"^  +  d^X^  +  ...  = 


0  (4.6) 


22'(/!)2 


V  =  M' 


(4.11) 
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2k 

dt  =  Sb,j,'  (4.12) 

i=k 

where  A  =  i-k. 

From  these  expressions  one  can  see  that  the  summation  notation  greatly  simplifies 
the  evaluation  of  eigenvalues. 

4.3.  Summary 

In  this  chapter,  a  technique  for  calculation  of  the  eigenv^ues  occurring  in  the  Graetz 
Problem  in  slip-flow  has  been  derived  by  constructing  a  matrix.  The  two  examples  show 
that  the  computation  of  the  bi,k  by  the  formulation  in  Eq.  (4.8)  is  much  more  effective  than 
when  found  by  the  recursion  relationship,  Eq.  (4.5).  With  the  formulation,  any  number  of 
eigenvalues  can  be  theoretically  determined.  The  next  chapter  will  deal  with  the  algo- 
nthms,  the  computational  program,  and  will  carry  out  the  calculation  of  the  eigenvalues. 


CHAPTER  5 


COMPUTATIONAL  RESULTS 


In  the  previous  chapter,  the  formulation  for  the  calculation  of  the  coefficients  occur¬ 
ring  in  the  eigenfunction  was  determined.  In  this  chapter,  we  will  develop  the  codes  for 
the  evaluation  of  eigenvalues  for  the  Graetz  Problem  in  slip-flow.  We  will  calculate  the 
eigenvalues  and  discuss  the  computational  results. 


5.1  IVeatment  of  Very  Large  Numbers  in  the  Computations 

By  using  Eq.  (4.8),  the  calculation  of  the  coefficient  I?,-/,  involves  the  constant 
which  can  become  a  very  large  number.  For  example,  if  i  =  30  (or  k  =  15),  then 
^  (50.0^  =  (7.75  X  10^^)(7.04  x  10^)  =  8.11  x  10^^ ,  which  will  overflow  on 

the  computer.  Therefore,  to  reduce  the  effect  of  veiy  large  numbers  in  the  computation, 
the  eigenfunction  must  be  treated  as  follows. 

(1)  We  can  combine  22'into  the  summation  in  Eq.  (4.8)  as 


^ijt  - 


22'(i!)2 


22^2"  +A-1)H 

s 


(-1)*  y  +  ^-1)^ 

sf=A  22 


=  y  +^-i)\  V  (5j_,  +/i-2)^ 

2"*-^(f!)2  s,=A  22  22 


So  for  i  -  30  (or  k  =  15) ,  the  largest  number  in  the  term  b2k,k  can  be  reduced  to  (30!)^  = 
7.04x10^. 
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(2)  Letting  A  _  X/g  and  d/^  _  g^^d^  and  allowing  the  magnification  coefficient,  g, 
be  any  given  number  greater  than  1 ,  we  have 

00 

G(X)  =  ^ ~  1  +  A'^^2  ••• 


=  .  +  + ...  = 


=  1  +  A'2j/  +X'%'  +X'%’  + 

00 

=  ^X'^dd  =  G{X')  =  0 


For  example,  when  we  let  g  be  10,  then  from  Eq.  (4.9)  and  Eq.  (4.8)  we  have 


=  ,^1020)^,. 

/=!  i=l 


1020)^,,  J  + 


di  =  =  102(2)^2  =  102(2)1:^,  =  2’l02(2)^., 

i=2  '  1=2 


=  102%,  +  102(2)Z,_+  102(2)a 


2(3)  2(3) 

ds'  =  =  102%  =  102(3)^^  =  ^102(3)^., 


=  102(3)Z,33  +  102(3)Z,43  +  102(3)^33  +  102(3)^^3 
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In  general  then, 


4'  = 


For  the  example  in  which  k  =  15,  we  have 

2(15)  2(15) 

d,,'  =  =  io2<'5)2:  fc  =  2'  I0WJ,„ 

i-15  i-is  ‘’I* 

Thus,  combining  lO^dS  into  the  tcnn  bio.,1 ,  the  largest  number  in  this  tern  reduces  to 

(i!)^/10^'^=  (30!)^/10^^=  7.04x  10^"^. 

(3)  The  magnitude  of  the  term  I  1  can  be  reduced  for  computational  purposes  by 
taking  the  logarithm  of  both  sides  and  later  reversing  this  computation  by  taking  the  expo¬ 
nential  function  of  both  sides.  That  is, 

=  -2i /ogjo2-2^X  logjok  ^^2; 

This  method  can  reduce  the  number  greatly.  For  example,  for  i  =  30,  then, 

^^^10^30.30  =  -2(30)  tog,o2-2^f  logjok 

=  -(60)  (0.301)-64.847  =  -82.907 

All  three  of  these  methods  were  used  in  the  computer  code  to  accommodate  the  in- 
herently  large  numbers. 

5.2  Flow  Chart  of  Computation 

A  block  diagram  for  computation  of  the  eigenvalues  is  shown  in  Figure  5.1.  The 
input  data  includes  the  number  of  ternis,  k.  the  magnification  coefficient,  g,  in  Eq.  (5.1), 
and  the  Knudsen  number,  Kn  oryS  in  Eq.  (4.11). 
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The  calculation  procedure  can  be  broken  down  into  three  steps: 

1)  calculation  of  coefficient  bij  (which  requires  a  summation  solver  for  Eq.  (4.8) 
for  bij  and  Eq.  (4.11)  for  bij  ’); 

2)  calculation  of  4  by  Eq.  (4.9)  and  4  ’  by  Eq.  (4.12);  and 

3)  calculation  of  eigenvalues  or  by  Eq.  (5.1). 

Based  on  the  flow  chart,  the  codes  for  computation  of  the  eigenvalues  have  been  de¬ 
veloped  and  are  listed  in  Appendix  A.  The  codes  need  as  input  the  value  of  the  Knudsen 
number,  Kn,  and  the  number  of  coefficients  4.  ork,  and  the  magnification  coefficient,  g. 
The  output  of  the  calculation  includes  the  coefficient  bij  ,  4,  and  eigenvalues, 

From  Examples  1  and  2  in  Chapter  4  for  the  formulation  of  Eq.  (4.8),  we  can  see  the 
feature  of  the  calculation  of  the  k-\h  term:  the  calculation  involves  all  the  previous  (k-l) 
terms,  that  is,  the  summations  of  parameters  sj,  S2. ...,  Sfe_],  and  further  expanded  summa¬ 
tions  of  si,  S2, ...,  Sk-i ,  for  the  present/:.  This  concept  can  be  illustrated  using  the  following 
triangles  in  Fig.  5.2. 


Fig.  5.2  Dlustration  of  the  growth  of  si,  S2, ... ,  Sk.i  as  functions  of  it 


5.3  Results 

5.3.1  Comparison  of  the  Firsf  Tph 
dir(Kn=0) 

5^3.1.1  Accuracy.  Because  the  evaluation  of  eigenvalues  involves  the  determination  of 
an  infinite  number  of  terms,  it  is  extremely  difficult  to  determine  directly  the  accuracy  of 
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the  approximate  numerical  eigenvalues.  All  we  have  to  do  is  to  find  the  consistencies 
existing  in  the  results  by  different  methods  instead  to  determine  the  relati\  e  accuracy. 

In  order  to  assess  the  accuracy  of  the  computation  of  Eq.  (4.8),  a  comparison  of  the 
first  ten  <4  using  Mathcad  5.0  was  made,  as  shown  in  Table  5. 1 .  In  the  table,  the  data  in  the 
second  column  are  the  computational  coefficients  4  (  with  g  =  10)  using  Eq.  (4.8);  the 
data  in  the  fourth  column  are  the  absolutely  exact  coefficients  4  calculated  by  expansions 
using  the  symbolic  processor  in  Mathcad  5.0;  the  third  column  gives  the  equivalent  deci¬ 
mal  values  of  the  fourth  column.  The  differences  between  Eq.  (4.8)  and  the  exact  values 
are  shown  in  Table  5.2. 


Table  5.1  Comparison  of  Coefficients  4  (^=10) 


k 

Eq.  (4.9) 

Simulation  by  Mathcad  Version  5.0 

Numeric 

equivalent 

Symbolic  solution 

0 

i.OOOO 

1.00000 

1 

I 

-18,7500 

-18.75000E-02 

-3/16 

2 

79.2101 

79.21 0069E-04 

73/9,216 

3 

-144.043 

-144.04297E-06 

-59/409.600 

4 

145.080 

145.07980E-08 

603.793/416,197,814,400 

5 

-92.6715 

-92.67 144E-10 

-555,379/59,929,893,273.600 

6 

40.8619 

40.861856E-12 

4.266,870,48 1/1 04,42 1 ,846,039,920,640,000 

7 

-13.1812 

-13.181156E-14 

-37.2 1 7,872, 1 47/282,356, 67 1 .69 1 ,945,4 1 0,560.000 

8 

3.24941 

3.24513 15E-1 6 

41,377.942.693.441/127,507,755,229,335,476.282.327.040.000 

9 

-0.646315 

-0.62971547E-18 

-9,28 1 ,940,782,645,85 1/1 4,739,896,504,5 1 1 ,1 81 ,058,237,005,824E+6 
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Table  5.2  Differences  of  4  Between  Two  Methods 


k 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Difference 

0 

1 

3.1E^9 

3.0E-11 

2.0E^12 

5.5E-15 

1.4E-17 

4.4E-19 

4.3E-19 

1.7E-20 

From  the  comparison,  we  can  see  that  the  coefficient  4  computed  by  Eq.  (4.8)  is 
always  a  little  larger  than  that  found  by  Mathcad,  but  the  differences  are  rather  small. 
Based  on  the  comparison,  it  is  evident  that  the  computational  results  have  enough  accura¬ 
cy  for  an  accurate  determination  of  the  eigenvalues. 

5,3t1,2  Efficiency.  In  order  to  assess  the  efficiency  of  the  computation  of  Eq.  (4.8),  we 
can  compare  the  computational  efficiencies  using  Mathcad  and  a  calculation  of  Eq.  (4.8) 
using  a  Fortran  code  for  it  =  9. 

Using  Mathcad  5.0,  three  steps  were  implemented  to  solve  for  the  eigenvalues: 

1)  Expand  a^,  a\,  a2, ...,  uig,  as  shown  in  Chapter  4; 

2)  Add  up  the  terms  in  the  expansions  with  respect  to  for  d],  X'^  for  d2, ... , 

X^^  for  dg;  and 

3)  Solve  the  eigenfunction  equation  for  the  eigenvalues 

1  -f  djX^ d2X‘^  +  d3X^  4- ...  d9X^^  =  0 

Procedure  1 )  took  about  ten  hours  and  2)  about  five  hours  to  accomplish.  For  large 
k,  for  example  k  =  25,  it  would  be  too  time-consuming  and  complex  to  use  those  methods. 

Using  the  computer  code  given  in  Appendix  A,  which  uses  Eq.  (4.8)  to  determine  all 
4  for  X  =  1  to  9,  the  process  took  approximately  one  minute.  Therefore,  the  latter  proce¬ 
dure  IS  an  effective  and  efficient  technique  for  the  calculation  of  eigenvalues. 
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5.3.2  Behavior  of  the  Eigenfiinrfinn 

In  practice,  the  number  of  terms  used  in  the  eigenfunction  expression  must  be  lim¬ 
ited.  In  order  to  see  how  the  eigenfunction  behaves  with  respect  to  the  number  of  coeffi¬ 
cient  terms,  a  number  of  calculations  have  been  carried  out  and  plotted. 

Figure  5.3  shows  the  behavior  of  the  eigenfunction  as  a  function  of  the  number  of 
coefficient  terms.  The  numbers  of  the  curves  indicate  the  number  of  terms  4  used  in  the 
computations.  At  least  six  terms  are  needed  to  obtain  the  two  lowest  eigenvalues,  shown 
as  the  curve  numbered  6. 


G(X’) 


Fig.S.3  Behavior  of  the  eigenfunction  as  the  number  of  coefficient  terms  increases 

Graetz  found  the  first  two  values,  X,  =2.7043  and  Xj  =  6.50  from  the  following  eigen- 
function: 

0  =  1  -  X2  0.1875  +  0.007921  -  X^  0.00014404  +  X^  145.92x10-8  + 

-  Xio  94.938x10-10-1-  ...  rc  .v 
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Comparing  the  coefficients  4  from  Eq.  (5.3)  with  those  given  in  Table  5. 1  we  can  see 
significant  differences  exist  from  the  5th  term  on,  revealing  why  Graetz’s  second  eigen¬ 
value  is  not  accurate. 

Obviously,  the  second  eigenvalue  determined  by  curve  6  is  less  than  the  ’’accurate 
value,”  which  is  slightly  greater  than  6.6  (see  Table  5.3),  shown  as  curves  7, 8,  and  9  etc., 
in  Figure  5.4. 

In  Figure  5.4,  we  can  see  that  the  true  second  eigenvalue  lies  between  those  obtained 
by  curves  7  and  8.  It  is  clear  that  the  second  eigenvalue  obtained  initially  by  curve  6  is 
rather  rough;  with  one  more  term,  the  value  is  larger  than  the  first  one  as  shown  by  curve  7, 
but  closer  to  the  accurate  value”;  the  second  eigenvalue  given  by  curve  8  is  a  little  less 
than  that  given  by  curve  7  but  is  even  closer  to  the  ’’accurate  value.”  The  same  holds  true 
for  the  case  with  nine  terms,  ten  terms,  and  so  on. 

From  this  discussion  we  can  conclude: 

(1)  that  the  eigenvalues  initially  obtained  with  the  minimum  number  of  coeffi¬ 
cients  is  always  rather  rough,  such  as  the  second  eigenvalue  obtained  by  curve 
6  and  the  third  eigenvalue  obtained  by  curve  9; 

(2)  that  the  next  to  the  last  available  eigenvalue  that  can  be  determined  for  a  certain 
number  of  coefficients  is  always  correct  and  sufficiently  accurate.  For  instance, 

for  curve  7,  the  second  eigenvalue  can  be  assumed  to  be  reasonably  accurate  and 
correct;  and 

(3)  that  the  eigenvalues  and  the  convergence  of  the  eigenfunction  are  sensitive  to  the 
accuracy  of  the  coefficients  4. 

Figure  5.5  shows  the  plot  of  the  eigenfunction  with  25  terms.  There  are  six  eigenva¬ 
lues  shown  in  the  plot.  From  the  above  discussion,  we  conclude  that  the  first  five  eigenva¬ 
lues  are  correct,  but  the  last  one,  or  the  sixth,  is  somewhat  inaccurate  due  to  the  truncation 


G(X’) 
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Gao 


Fig.  5.5  Plot  of  eigenfunction  with  25  coefficient  terms 


of  the  eigenfunction  expression.  The  plot  also  shows  that  the  eigenfunction  is  oscillating 
with  decreasing  magnitude. 

£3.3  Comparison  With  Previmicly  Known 
Eigenvalues  With  y/i  =  0 

Table  5.3  shows  the  comparison  with  previously  known  eigenvalues  for  the  classical 
’’Graetz  Problem”  (^  =  1  or  /ifn  =  0)  as  presented  by  Sellars  et  al.  ( 1 956).  As  seen  in  Table 
5.3,  the  first  four  eigenvalues  are  in  excellent  agreement.  Therefore,  we  can  apply  this 
technique  to  the  Graetz  Problem  in  slip-flow  for  evaluation  of  the  eigenvalues. 

5.3.4  Eigenvalues  for  ATw  n 

Table  5.4  shows  the  first  five  eigenvalues  for  slip-flow  with  different  Knudsen  num¬ 
bers,  Kn. 


Table  5.3  Comparison  with  Previously  Known  Eigenvalues 
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G(X’) 


Figure  5.6  Plots  of  the  eigenfunction  as  a  function  of  Knudsen  number 


Figure  5.6  shows  the  behavior  of  the  eigenfunction  for  various  Kn  under  slip-flow  condi¬ 
tions.  It  shows  that  the  eigenvalues  decrease  as  Kn  increases.  For  Kn  >  0,  the  plots  appear 
unstable  after  the  fifth  root  so  that  only  the  first  four  values  are  reliable.  The  possible 
cause  for  this  instability  is  that  the  truncation  errors  are  magnified  by  the  factor  (l-h4A'/i)' 

on  m  the  modified  matrix  B.  The  coefficients  4  ’  for  Kn  from  0.00  to  0. 12  are  shown 
in  Appendix  B. 

5,3.5  Influence  of  Kn  on  the  Nnsspit  Niimhpr 


£3.5.1  Local  heat  transfer  coefficient.  Using  Eq.  (3.29)  and  Eq.  (3.20),  the  local  heat 
transfer  coefficient  Nux  has  been  calculated. 
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Figure  5.7  shows  the  local  Nux  value  as  a  function  of  x*/Gz  for  Kn  =  0.02  and  with  the 
number  of  eigenvalues  as  a  parameter.  The  value  of  the  local  Nusselt  number  converges 
dramatically  with  the  increase  in  the  number  of  eigenvalues  in  the  computation.  When 
x*/Gz  is  0.02,  the  error  in  Nu^  is  0.7  percent  when  two  eigenvalues  are  used  and  compar¬ 
ing  to  the  straight  line  ( using  one  eigenvalue ),  the  error  is  1 4  percent.  It  can  be  concluded 
that  the  results  using  four  eigenvalues  are  sufficiently  accurate  for  x*/Gz  >  0.02.  When 
x*/Gz  is  greater  than  0.05,  the  error  is  at  most  1 .3  percent — that  is,  all  three  plots  become 
nearly  flat,  indicating  a  thermally  fully-developed  condition. 


Figure  5.8  shows  the  local  Nusselt  numbers  as  a  function  of  Kn.  It  is  obvious  that  Kn 
has  an  influence  on  the  Nusselt  number.  All  the  plots  in  Fig.  5.8  show  that  the  Nusselt 
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number  increases  as  Kn  increases,  and  that  this  effect  is  magnified  near  the  entrance. 
When  x*/Gz  is  greater  than  0.05,  all  the  plots  become  nearly  flat,  indicating  a  thermally 
fully-developed  condition. 


Nu(x) 


Fig.  5.8  The  local  Nusselt  numbers  as  functions  of  x*/Gz  and  Kn  for  n  =  4 

Table  5.5  shows  the  values  of  the  Nusselt  number  for  fully-developed  conditions 
(where  x*/Gz  >  0.05)  for  different  Kn  and  the  ratios  of  these  values  to  those  with  Kn  =  0. 
This  ratio  increases  with  an  increase  of  Kn.  The  data  show  that  when  Ati  is  0.0 1 ,  the  value 
of  the  Nusselt  number  increases  about  3  percent,  and  when  Kn  is  0.02,  the  increase  is 
greater  than  5  percent.  Thus,  we  can  conclude  that  when  Kn  is  greater  than  0.0 1 ,  the  effect 
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Fig,  5.9  Fully  developed  Nu  as  a  function  of  Kn 

5.3.5.2  Overall  heat  traasfer  coefficient.  Using  Eqs.  (3.37),  (3.34),  and  (3.39),  the  over¬ 
all  heat  transfer  coefficient  Nu  can  be  calculated.  Fig.  5.10  shows  the  plots  of  i^,  Nu^, 
and  (NiT-  Nuoo),  or  A  as  functions  of  x*/Gz. 

From  Fig.  5. 1 0,  we  can  see  that  the  overall  heat  transfer  coefficient  Nu  is  greater  than 
the  local  heat  transfer  coefficient  Nux:  that  the  entrance  has  a  greater  effect  on  Nu  than 
on  Nux:  and  that  when  x*/Gz  is  greater  than  0.05,  Nu^  becomes  nearly  flat,  indicating  a 
thermally  fully-developed  condition.  However,  Nu  does  not  becomes  flat,  so  that  it  can 
not  be  considered  as  a  fully-developed  condition  until  x*/Gz  greater  than  0.20. 

From  the  above  discussion,  we  can  see  that: 

( 1 )  slip-flow  has  a  positive  influence  on  the  heat  transfer  coefficient  and  can  enhance 
the  heat  transfer  efficiency; 
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Fig.  5.10  Plots  of  Nu,  Nux  and  A  as  functions  of  dimensionless  axial  location 


(2)  the  influence  depends  on  the  Knudsen  number  and  increases  as  Kn  increases; 

(3)  when  Kn  is  equal  or  greater  than  0.02,  the  increase  in  the  fully-developed  Nu  is 
greater  than  5  percent  so  that  the  effect  of  slip— flow  should  be  taken  into  consider¬ 
ation  in  the  computations  of  the  heat  transfer  coefficient;  and 

(4)  that  the  influence  of  Kn  on  Nuoo  will  decrease  as  Kn  increases. 

5.4  Simplified  Eigenvalue  Relationships 
with/i:/i 

Figure  5.11  shows  that  are  functions  of  Kn.  For  practical  purposes,  a  simplified 
expression  for  calculation  of  the  eigenvalues  is  needed.  By  using  a  least-squares  curve 
fit  program,  the  following  exponential  expressions  as  given  in  Table  5.4  were  found  to 
yield  the  best  fit.  The  general  form  may  be  expressed  as 
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5.5  Summary 

In  this  chapter,  the  codes  were  developed  for  the  evaluation  of  eigenvalues.  The  first 
six  eigenvalues  were  found,  and  the  behavior  of  the  eigenfunction  was  plotted.  From  the 
comparison  and  discussion,  it  is  evident  that  the  new  technique  using  Eq.  (4.8)  for  evalua¬ 
tion  of  the  eigenvalues  of  the  Graetz  Problem  in  slip-flow  is  computionally  effective  and 
efficient;  the  Nusselt  number  increases  as  the  Knudsen  number  increases;  and  the  simpli¬ 
fied  relationship  of  the  eigenvalues  as  a  function  of  the  Knudsen  number  is  reliable  and 
convenient  for  calculation  purposes. 


CHAPTER  6 


CONCLUSIONS  AND  FURTHER  RESEARCH 

6.1  Conclusions 

In  the  previous  chapters,  the  mathematical  models  of  velocity  distribution  and  tem¬ 
perature  distribution  were  established,  and  the  expression  for  the  series  solution  shows 
the  importance  of  the  eigenvalues.  Since  those  eigenvalues  were  extremely  difficult  to 
evaluate  direcUy  from  the  original  expansion,  a  formulation  was  derived  based  on  a  spe- 
ciaUy  constructed  matrix  of  coefficients  bi,j.  The  formulation  can  be  used  to  find  any  giv¬ 
en  bij  using  only  the  indices  i  and/  TTiis  fact  makes  it  possible  to  evaluate  the  eigenvalues 
by  computer.  The  computer  codes  were  developed  and  some  results  were  obtained.  From 

the  discussions  and  analysis  of  the  computational  results,  the  following  conclusions  can 
be  drawn: 

1 .  The  technique  for  evaluation  of  the  eigenvalues  of  the  Graetz  Problem  in  slip-flow 
is  computational  effective; 

2.  The  Nusselt  number  increases  as  Kn  increases,  or  the  heat  transfer  is  enhanced 
under  slip-flow  conditions; 

3.  When  Kn  is  equal  to  or  greater  than  0.02,  the  increase  in  fully  developed  Nusselt 
number  is  greater  than  5  percent  so  that  the  effect  of  slip  flow  conditions  should 
be  taken  into  consideration  in  the  computations  of  the  heat  transfer  coeffient;  and 

4.  The  simplified  relationship  between  the  eigenvalues  and  the  Knudsen  number  is 
reliable  and  convenient  for  calculation  purposes. 
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6.2  Further  Research 

The  evaluation  of  the  eigenvalues  is  veiy  important  for  the  solution  of  the  Graetz 
Problem  in  slip-flow.  Although  the  technique  is  effective  for  n  <  5 ,  it  is  extremely  time- 
consuming,  and  the  computational  eiror  is  a  problem  for  large  1.  Based  on  this  work, 
the  following  future  research  is  suggested: 

1.  Develop  computer  codes  for  calculating  the  heat  transfer  coefficient  and  Nusselt 

number,  which  involves  the  computation  of  Q  and  differentiation  of  eigenfunction 
G„  atr*=  1; 

2.  Obtain  a  simplified  relationship  between  the  overall  Nusselt  number  and  the 
Knudsen  number; 

3.  Develop  a  more  effective  technique  to  deal  with  very  large  numbers  in  computa¬ 
tions; 

4.  Improve  the  codes  to  reduce  the  computing  time;  and/or 

5.  Develop  a  more  effective  technique  to  reduce  the  computing  time. 
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♦**♦**+♦*♦******„**+****», 

*  THIS  PROGRAM  FOR  CALCULATION  OF  COEFRCIENT  a;  j  and  d^ 

*  a:  aij 

*  at  initial  of  aij  * 

*  d:  dk  » 

*  g:  magnification  coefficient  ♦ 

*  kn:  Knudsen  number  * 

*  bb:  sqrt(L+4.*kn)  « 

*  s:  A  =  i-j  * 


real  a(  1 000, 1 000), d(0: 1 000),r,u,at,sum,g,kn,bb 
integer  s 

C  INPUT  Kn,g 

read  *,kn,  g 
bb=sqrt(  1  .+4.  *kn) 

C  CALCULATING  aij 
do  10i=l,26 

a(i,i)=(-l)**i*(((g/2)**(i)/(funtl(i)))**2)*bb**(2*i) 

10  continue 

C  CALCULATING  aij 
do  50j=I,25 
at=l 

print  *,’*** ’,j 
do40i=j+l,2*j 
s=i-j 

summ(i,j,s,at) 

a(i,j)=at*(-l)**i*(((b/2)**(j)/((funtl(i)/10.**(s/2.)/(2**s))))**2)*bb**(2*i-4*s) 

C  OUTPUT  aij 

write(*,I10)i,j,a(i,j) 

1 1 0  fonnat(L5,5x,i5,5x,e  1 0.4) 

at=l 

40  continue 

*  print*,a(i,j) 

50  continue 

C  INrnATlNG  dk  zreo 
do  55j=0,25 
d0)=0 

55  continue 


C  CALCULATING  dk  BY  SUMMATION  aicj 
d{0)=l 
do60j=l,25 
do  70  i=j,2*j 

d(j)=d{j)+a0d) 

70  continue 

60  continue 

C  OUTPUT  k,Kn,dk 

open(3) 

write(3,*)’k  =  25  Kn  =  ’,kn 
do  100  i=0.25 
print  ♦,d(i) 
write(3,*)d(i) 

100  continue 

C  CALCULATING  THE  VALUES  OF  EIGENFUNCTION  AT  EACH  POINT 
du=0.0001 
open(3) 

do  200  i=  1,30000 
u=du*(i-l) 

write(3,*)u,f(u,25,d),df(u,25,d) 

200  continue 

call  newton(d,25) 
end 

C  FUNCTION  OF  n! 

function  funtl(n) 
fiintl=l 
do  10  i=l,n 
funtl=funtl’'‘i 
10  continue 

return 
end 

*(1)  function  OF  THE  LOWEST  SUMMATION 

function  suml(jj) 
real  suml 

suml=0 
do  10i=l,jj 
suml=suml+i**2/10. 

10  continue 

suinl=sunil/4. 

return 

end 


*(2)  function  OF  THE  2ND  LOWEST  SUMMATION 

function  suin2(jj) 
real  suin2 
integer  ss 
sum2=0 
do  10  ss=2,j|j 

sum2=sum2+suml  (ss-1  )*(ss+ 1  )**2/1 0. 

10  continue 

sum2=suin2/4. 

return 

end 

♦(3) 

function  sum3(jj) 
real  sum3 
integer  ss 

sum3=0 
do  10  ss=3,j[j 

sum3=sum3+(ss+2)**2*sum2(ss-l)/10. 

10  continue 

sum3=sum3/4. 

return 

end 

*(4) 

function  sum4(jj) 
real  sum4 
integer  ss 

suin4=0 
do  10  ss=4,j[j 

sum4=suin4+(ss+3)**2*sum3(ss-l)/10. 

10  continue 

sum4=sum4/4. 

return 

end 

*(5) 

function  sum5(jj) 
real  sum5 
integer  ss 

sum5=0 
do  10  ss=5,ij 

suni5=suin5+(ss+4)**2*suin4(ss-l)/10. 

10  continue 


siim5=sum5/4. 

return 

end 

*(6) 

function  suin6(jj) 
real  suin6 
integer  ss 
sum  6=0 
do  10  ss=6,jy 

sum6=sum6+sum5(ss-l)*{ss+5)**2/10. 

10  continue 

sum6=sum6/4. 

return 

end 

*(7) 

function  sum7(jj) 
real  sum7 
integer  ss 

sum7=0 
do  10ss=7,ij 

sum7=sum7+(ss+6)*  *2*sum6(ss-l  )/l  0. 

10  continue 

sum7=sum7/4. 

return 

end 

*(8) 

function  sum8(jj) 
real  sum8 
integer  ss 

sum8=0 
do  10ss=8,j|j 

sum8=sum8+(ss+7)**2*sum7(ss-l)/10. 

10  continue 

suin4=suin4/4. 

return 

end 

*(9) 

function  sum9(jj) 
real  sum9 
integer  ss 
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sum  9=0 
do  10  ss=9,jy 

$iim9=sum9+(ss+8)*  *2*siim8  (ss- 1 )/ 1 0. 

10  continue 

sum9=sum9/4. 

return 

end 

♦  ♦  ♦  ♦  *  3|t  :<t  ♦  ♦  +  3fe  *  *  ♦  *  3*^  j|t  3(c  *  ^  *  jjc  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

*(10) 

function  sumlO(jj) 
real  sum  10 
integer  ss 

sum  10=0 
do  10  ss=10,jj 

sum  1 0=sum  1 0+(ss+9)  **2*sum9(ss--l )/ 1 0. 

10  continue 

sum  10=sum  10/4. 

return 

end 

*(11) 

function  sumll(jj) 
real  sumll 
integer  ss 

sumll=0 
do  10  ss=ll,jj 

sum  1 1  =sum  1 1  +(ss+ 1 0)*  *2*sum  1 0(ss-l  )/l  0. 

10  continue 

sumIl=sumll/4. 

return 

end 

♦(12) 

function  suml2(jj) 
real  sum  12 
integer  ss 

sum  12=0 
do  10  ss=12,jj 

sum  12=sum  1 2+(ss+ 1 1  )**2»sum  1 1  (ss-1  )/10. 

10  continue 

sum  1 2=sum  1 2/4. 

return 

end 

*(13) 


function  suml3(jj) 
realsumlS 
integer  ss 

sum  13=0 
do  10  ss=13,ij 

sum  1 3=s  urn  1 3+(ss+ 1 2)  *  *2  ♦sum  1 2(ss- 1  )/l  0. 
10  continue 

sum  13=sum  13/4. 

return 

end 

*(I4) 

function  suml4(jj) 
real  sum  14 
integer  ss 

suml4=0 
do  10  ss=14,ij 

suml4=suml4H-(ss+13)**2*suml3(ss-l)/10. 

10  continue 

sum  14=sum  14/4. 

return 

end 

function  suml5(jj) 
real  suml5 
integer  ss 
sum  15=0 
do  10  ss=15,jy 

sum  1 5=sum  1 5+(ss+ 1 4)**2*sum  14{ss-l  )/l  0. 
10  continue 

sum  15=sum  15/4. 

return 

end 

*(16) 

function  sumI6(jj) 
real  suml6 
integer  ss 

suml6=0 
do  10ss=16,ij 

suml6=suml64-(ss+15)**2*suml5(ss-l)/10. 

10  continue 

sum  16=sum  16/4. 
return 


*(17) 

function  suml7(jj) 
real  suinl7 
integer  ss 

sum  17=0 
do  10  ss=17,ij 

sum  1 7=sum  1 7+(ss4- 1 6)  *  *2*sum  1 6(ss~l  )/l  0. 
10  continue 

sum  17=sum  17/4. 

return 

end 

♦  3(c  3(e Jjt  + :ft  ♦♦  **  :(t  :ic  5jc  **  :4c  3)£  :jt  :)e  *  4;  *♦  5(s  *  3^  ^  :(t  3je  ******  + 

*(18) 

function  suml8(jj) 
real  suml$ 
integer  ss 

sum  18=0 
do  10  ss=18,j[j 

suml8=suml8+(ss+17)**2*suml7(ss-l)/10. 
10  continue 

sum  18=sum  18/4. 

return 

end 

*(19) 

function  suml9(jj) 
real  sum  19 
integer  ss 

sum  19=0 
do  10ss=19,jj 

suml9=sum  1 9+(ss+ 1 8)**2*sum  1 8(ss-l  )/l  0. 

10  continue 

suml9=suml9/4. 

return 

end 

*(20) 

function  sum20(jj) 
real  sum20 
integer  ss 

sum20=0 
do  10  ss=20,ij 

sum20=sum2Of(ss+ 1 9)**2*suml  9(ss-l)/10. 


c  print*, ’suin20=’,suin20 

10  continue 

sum20=suin20/4. 

return 

end 

*(21) 

function  sum21(jj) 
real  suin21 
integer  ss 

sum21=0 
do  10ss=21,ij 

sum21=suni21+(ss+20)**2*suin20(ss-l)/10. 
c  print*,’sum21=’  ,sum2I 

10  continue 

sum21=suni21/4. 

return 

end 

*(22) 

function  sum22(jj) 
real  suin22 
integer  ss 

sum22=0 
do  10  ss=22,ij 

sum22=sum22+(ss+2 1  )**2*suin2 1  (ss-1  )/10. 
c  print*, ’suin22’,suni22 

10  continue 

sum22=suin22/4. 

return 

end 

♦  ♦  :(t  :)e  :jc  j(c  :(e  ♦  3<e  :|e  5je  5|c  3ft 3jc :({ jJ:  5|c  :^e  ;je  jjc :({ 3(c  3^t  *  j(j  5|c  *  jfe  *  :(c  + 4: 

*(23) 

function  suin23(jj) 
real  suin23 
integer  ss 

sum23=0 
do  10  ss=23,ij 

sum23=sum23+(ss+22)**2*sum22(ss-l  )/l  0. 

10  continue 

sum23=sum23/4. 

return 

end 

*(24) 


function  sum24(ij) 
real  suiii24 
integer  ss 


suni24=0 
do  10  ss=24,j[j 

suin24=sum24+(ss+23)*  *2*sum23(ss-l  )/l  0. 

10  continue 

sum24=sum24/4. 

return 

end 

*(25) 

function  sum25(ij) 
real  suin25 
integer  ss 
sum25=0 
do  10  ss=25,ij 

sum25=sum25+(ss+24)**2*suin24(ss-l)/10. 

10  continue 

suin25=sum25/4. 

return 

end 

subroutine  newton(para,k) 
real  para{0:1000),k 
read  *,ul,pr 

10  u2=ul-f(ul,k.para)/df(ul4c,para) 

error=abs((u2-ul)/u2) 
if(erTor.gtpr)  then 
ul=u2 

print  *,01, error 
goto  10 
endif 
u2=b*u2 

print  *,u2,f(u2/2,k,para) 

return 

end 

♦  ♦*♦♦******:(. 

function  f(udc,para) 
real  parafO:  1000) 
f=0. 

dol0i=0dc 

f=f+para(i)*u**(2*i) 

10  continue 

return 


end 

function  df(u4c,para) 
real  para(0:l(X)0) 
df=0. 

do20i=14c+l 

df=df+para(i)*u**(2*i-l)*2*i 
20  continue 

return 

C  FUNCTION  OF  SUMMATION 
function  summ(i,j,s,at) 
real  at 

if  (s  .eq.  24)  then 

at=sum24(j) 
goto  30 

endif 

if  (s  .eq.  25)  then 

at=sum25(j) 
goto  30 

endif  if  (s  .eq.  23  )  then 
at=sum23(j) 
goto  30 

endif 

if  (s  .eq.  22)  then 

at=sum22(j) 
goto  30 

endif 

if  (s  .eq.  21)  then 

at=sum21(j) 
goto  30 

endif 

if  (s  .eq.  20)  then 

at=sum20(j) 
goto  30 

endif 

if  (s  .eq.  19)  then 

at=suinl9(i) 
goto  30 

endif  if  (s  .eq.  18  )  then 
at=suml8(j) 
goto  30 

endif 

if(s  .eq.  17)  then 

at=suml7(j) 
goto  30 
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endif 

if  (s  .eq.  16)  then 

at=suml6(j) 
goto  30 

endif 

if  (s  .eq.  14)  then 

at=suinl4(j) 
goto  30 

endif 

if  (s  .eq.  15)  then 

at=suml5(J) 
goto  30 

endif  if  (s  .eq.  13  )  then 
at=suml3(j) 
goto  30 

endif 

if  (s  .eq.  12)  then 

at=suml2(j) 
goto  30 

endif 

if  (s  .eq.  11)  then 

at=sumll(j) 
goto  30 

endif 

if  (s  .eq.  10)  then 

at=suml0(j) 
goto  30 

endif 

if  (s  .eq.  9)  then 

at=sum9(j) 
goto  30 

endif  if  (s  .eq.  8  )  then 
at=sum8(j) 
goto  30 

endif 

if(s  .eq.  7)  then 

at=sum7(j) 
goto  30 

endif 

if  (s  .eq.  6)  then 

at=sum6(j) 
goto  30 

endif 

if  (s  .eq.  4)  then 

at=sum4(j) 
goto  30 
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endif 

if  (s  .eq.  5)  then 

at=sum5(j) 
goto  30 

endif  if  (s  .eq.  3 )  then 
at=sum3(j) 
goto  30 

endif 

if(s  .eq.  2)  then 

at=sum2(j) 
goto  30 

else 

at=suml(j) 

endif 

30  returen 

end 
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THE  COEFFICIENT  4’  OF  EIGENFUNCTION  FOR  DIFFERENT  Kn 


k 

0.00 

0.02 

0.04 

0.06 

0.08 

0 

1.000000 

1.000000 

1.000000 

1.000000 

1.000000 

1 

-18.750000 

-20.750000 

-22.750000 

-24.750000 

-26.750000 

'2 

79.210100 

98.265600 

119.321000 

142.377000 

167.432000 

3 

-144.043000 

-201.296000 

-271.689000 

-356.553000 

^57.224000 

4 

145.080000 

229.024000 

344.257000 

497.614000 

696.594000 

5 

-92.671500 

-165.560000 

-277.502000 

-442.200000 

-676.202000 

6 

40.861900 

82.728000 

54.757000 

272.024000 

454.598000 

7 

-13.181200 

-30.273700 

-63.245900 

-122.683000 

-224.130000 

8 

3.249410 

8.466540 

19.757500 

42.302900 

84.501400 

9 

-0.646315 

-1.883630 

-4.887290 

-11.531000 

-25.168400 

10 

0.124692 

0.363041 

1.005790 

2.575100 

6.103080 

11 

-3.58121E-02 

-7.81416EM)2 

-0.196144 

-0.505606 

-1.259180 

12 

1.48968EM)2 

2.45660E-02 

4.68932E-02 

1.02985E-01 

0.243432 

13 

-5.81835E-03 

-9.07529E-03 

-1.47620E-02 

-2.62552E-02 

-5.21826E-02 

14 

1.82761E-03 

2.97430E-03 

4.78575E-03 

7.83189E-03 

1.34943E-02 

15 

-4.58083E-04 

-7.98782E-04 

-1.34809E-03 

-2.23449E-03 

-3.71120E-03 

16 

9.33666E-05 

1.75660E-04 

3.16460E-04 

5.51241E-04 

9.39192E-04 

17 

-1.57875E-05 

-3.21078&-05 

-6.20691E-05 

-1.15009E-04 

-2.05969E-04 

18 

2.25303E-06 

4.95574E-06 

1.02964E-05 

2.03759E^5 

3.869 13E-05 

19 

-2.75302E-07 

-6.55009E-07 

-1.46339E-06 

-3.09708E-06 

-6.255 17E-06 

20 

2.91544E-08 

7.50298E-08 

1.80276E-07 

4.08225E-07 

8.77967E-07 

21 

-2.70262E-09 

-7.52384E-09 

-1.94436E-08 

-4.71180E-08 

-1.07959E-07 

22 

222204E-10 

6.68238E-10 

1.85579E-09 

4.81032E-09 

1.17395E-08 

23 

-1.57586E-11 

-5.16599E-11 

-1.55104E-10 

-4.31820E-10 

-1.12584E-09 

24 

1.17557E-12 

4.049 16E-12 

1.28004E-11 

3.75388E^11 

1.03039E-10 

25 

-3.96765E-14 

-1.63832E-I3 

-6.01195E-13 

-2.00029E-12 

-6.12443E-12 
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THECOEFHCffiNTdk’  OF  EIGENFUNCTION  FOR  DIFFERENT  A:«  (continued) 


k  0.005 


0  1.000000000000 

1  -19.249999755282 

2  83.786455012016 

3  -157.19747929320 

4  163.48101833927 

5  -107.88202195819 

6  49.163092053819 

7  -16.395674045065 

8  4.1785060393633 

9  -0.85474908176110 

10  0.16295069004395 

11  -4.267884996578  lD-02 

12  1.66960731 644660-02 

13  -6.4933580354236D-03 

14  2.068888 1265729D-03 

15  -5.2835960492763D-O4 

16  1.09849683 15700D-04 

17  -1.8952355167967D-05 

18  2.7599 126784909D-06 

19  -3.4412736954670D-O7 

20  3.7186482821050I>-08 

21  -3.5185593897498D-09 

22  2.941 1030829215D-10 

23  -2.1886373955472D-11 

24  1.4598270425440D-12 

25  -8.7801704962411D-14 


0.10 


1.000000000000 
-28.750000000000 
194.48784722222 
-575.03255208333 
949.36357393887 
-999.17322697416 
728.53075977812 
-389.65391363012 
159.38657434817 
-51.472458492094 
13.456265606584 
-2.9059588237375 
0.52713284089290 
-8.1451719510272D-02 
1. 0849459 119075EM)2 
-1 .2586480084747D-03 
1.2831332514336D-04 
-1 . 1585667342762D-05 
9.3297867024356D-07 
-6.7424213213213rM)8 
4.397 125499401 7D-09 
-2.6008164814639D-10 
1.4015685304521D-11 
-6.91 0073243923  lD-1 3 
3.1286794400673D-14 
-1.3054533728857D-15 


0.12 


1.000000000000 

-30.750000000000 

223.54340277778 

-711.31380208333 

1264.7562007813 

-1434.1592495969 

1126.9391090331 

-649.68941116624 

286.48984269996 

-99.748319034885 

28.116494799992 

-6.5472351852919 

1.2806819641796 

-0.21339847213796 

3.06537181306020-02 

-3.83508782811710-03 

4.21646845579410-04 

-4.10594285963680-05 

3.56603668730940-06 

-2.77944923233250-07 

1.95499746628130-08 

-1.24717259466340-09 

7.24895093680860-11 

-3.85470738694550-12 

1.88244080230260-13 

-8.47181196833900-15 
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I.  bij  for  A'/z  =  0  (  P  =  1  ) 


j  I  i=j,j+l.j+2 . 2*j-l,2*j 

1  -.2500E+02  0.6250E+01 

2  0.1562E+03  -.8681E+02  0.9766E+01 

3  -.4340E+03  0.3798E+03  -.9657E+02  0.6782E+01 

4  0.6782E+03  -.8138E+03  0.3308E+03  -.5273E+02  0.2649E+01 

5  -.6782E+03  0.1036E+04  -.5840E+03  0.1498E-(-03  -.1709E+02  0.6623E+00 

6  0.4710E+03  -.8746E+03  0.6328E+03  -.2265E+03  0.4173E+02  -.3668E+01  0.1150E+00 

7  -.2403E-I-03  0.5256E+03  -.4653E+03  0.2148E+03  -.5541E+02  0.7881E+01  -.5597E+00 
0.1467E-01 

8  0.9386E+02  -.2364E+03  0.2477E+03  -.1405E+03  0.4692E+02  -.9363E-I-01  0.1077E+01 
-.6381E-01  0.5729E-O2 

9  -.2897E+02  0.8256E+02  -.lOOOE+03  0.6740E+02  -.2773E+02  0.7179E+01  -.1160E+01 
0.1114E+00  -.2258E-01  0.4420E-03 

10  0.7242E+01  -.2304E+02  0.3170E+02  -.2477E+02  0.1213E+02  -.3878E+01  0.8143E+00 
-.llOOE+00  0.3612E-01  -.1594E-02  0.2763E-04 

11  -.1496E+01  0.5258E+01  -.8098E+01  0.7202E+01  -.4100E+01  0.1564E+01  -.4059E+00 
0.7132E-01  -.3297E-01  0.2357E-03  -.9193E-04  0.1427E-05 

12  0.2598E+00  -.9992E+00  0.1704E+0I  -.1700E+01  0.1104E+01 -.4899E+00  0.1520E+00 
-.3305E-01  0.1989E-O1  -.2002E-02  0.1265E-03  -.4412E-05  0.6194E-07 

13  -.3843E-01  0.1606E+00  -.3003E+00  0.3324E+00  -.2425E+00  0.1230E+00  -.4445E-01 
0.1156E-01  -.8622E-02  0.1131E-02  -.1006E-03  0.5681E-05  -.1790E-06  0.2291E-08 

14  0.4902E-02  -.2211E-01  0.4501E-01  -.5473E-01  0.4436E-01  -.2531E-01  0.1046E-01 
-.3176E-02  0.2835E-02  -.4609E-03  0.5345E-04  -.4252E-05  0.2167E-06  -  6216E-08 
0.7304E-10 

15  -.5446E-03  0.2638E-O2  -.5806E-02  0.7696E-02  -.6862E-02  0.4355E-02  -.2028E-02 
0.7045E-03  -.7343E-03  0.1430E-03  -.2057E-04  0.2134E-05  -.1532E-06  0  7109E-08 
-.1869E-09  0.2029E^11 

16  0.5319E-04  -.2753E-03  0.6516E-03  -.9350E-03  0.9099E-03  -.6359E-03  0.3297E-03 
-.1291E-03  0.1542E-O3  -.3509E-04  0.6051E-05  -.7788E-06  0.7293E-07  -  4767E-08 
0.2027E-09  -.4916E-11  0.4953E-13 

17  -.4601E-05  0.2535E-04  -.6419E-04  0.9916E-04  -.1046E-03  0.7984E-04  -.4563E-04 
0.1991E-04  -.2683E-04  0.7000E-05  -.1411E-05  0.2178E-06  -.2532E-07  0  2158E-08 
-.1293E-09  0.5068E^11  -.1140E-12  0.1071E-14 

18  0.3550E-06  -.2074E-05  0.5596E-05  -.9260E-05  0.1052E-04  -.8716E-05  0.5446E-05 
-.2622E-05  0.3938E-05  -.1160E-05  0.2682E-06  -.4842E-07  0.6753E^8  -.7148E-09 
0.5585E-10  -.3084E-1I  .1121E-12  -.2348E-14  0.2066E-16 


19  -.2459E-07  0.1518E-06  -.4348E-06  0.7674E-06  -.9350E-06 
0.835  lE-06  -.5665E-06  0.2984E-06  -.4947E-06  0.1626E-06 
-.4246E-07  0.8795E-08  -.1435E-08  0.1823E-09  -.1769E-10 
0.1274E-11  -.6523E-13  0.2207E^14  -.4326E^16  0.3578E-18 

20  0.1537E-08  -.lOOOE-07  0.3030E-07  -.5681E-07  0.7389E-07 
-.7080E-07  0.5183E-O7  -.2966E-07  0.5383E-07  -.1954E-07 
0.5695E-08  -.1333E-08  0.2498E-09  -.3713E-10  0.4325E-11 
-.3871E-12  0.2584E-13  -.1231E-14  0.3896E-16  -.7169E-18 
0.5590E^20 

21  -.8711E-10  0.5959E-O9  -.1905E-08  0.3781E-08  -.5228E-08 
0.5351E-O8  -.4206E-08  0.2599E-08  -.5127E-08  0.2038E-08 
-.6564E-09  0.1716E-09  -.3634E-10  0.6204E-11  -.8461E-12 
0.9090E-13  -.7540E-14  0.4687E-I5  -.2088E-16  0.6200E-18 
-.1075E-19  0.7922E-22 

22  0.4499E-1I  -.3228E-10  0.1085E-09  -.2274E-09  0.3332E-09 
-.3628E-09  0.3047E-09  -.2023E-O9  0.431  IE-09  -.1864E-09 
0.6578E-10  -.1901E-10  0.4497E-11  -.8683E-12  0.1360E-12 
-.1711E-13  0.1704E^14  -.1316E^15  0.7650E-17  -.3198E-18 
0.8939E-20  -.1464E-21  0.1023E-23 

23  -.2126E-12  0.1596E-11  -.5632E-11  0.1242E-10  -.1922E-10 
0.22I9E^10  -.1984E-10  0.1408E-10  -.3226E-10  0.1508E-10 
-.5789E-11  0.I834E-11  -.4797E-12  0.1035E-12  -.1834E-13 
0.2651E-14  -.3093E-15  0.2869E-16  -.2071E-17  0.1130E-18 
-.4444E-20  0.1173E-21  -.1820E-23  0.1209E-25 

24  0.9229E-14  -.7235E-13  0.2673E-12  -.6192E-12  0.1009E^11 
-.123IE-11  0.1168E-11  -.8826E-12  0.2163E-11  -.1087E-11 
0.4512E-12  -.1556E^12  0.4463E-13  -.1065E-13  0.2109E-14 
-.3450E-15  0.4625E-16  -.5026E-17  0.4359E-18  -.2954E-19 
0.1516E-20  -.5632E-22  0.1408E-23  -.2073E-25  0.1312E-27 

25  -.3692^15  0.3017^14  -.I165E-13  0.2827E-13  -.4841E-13 
0.6223E-I3  -.6241E-13  0.5006E-13  -.1307E-12  0.7028E-13 
-.3139E-13  0.1171E-13  -.3656E-14  0.9570E-15  -.2097E-15 
0.3835E-16  -.5819E-17  0.7267E-18  -.7385E-19  0.601  lE-20 
-.3834E-21  0.1858E-22  -.6533E-24  0.1550E-25  -.2172E-27 
0.13I2E-29 


APPENDIX  D 


THE  INTEGRAL  OF  EQ.  (3.18) 


The  constant  Q  in  the  solution 


e  {r*,x*)  =  ^  C„G„(r*)  exp  [  -  J:*  (1  +  8  Afn) 

n  =  i  Gz 

must  satisfy  the  following  condition  for  all  r* 

00 

1  “  ^  CfiGfi 


(D.l) 


(D.2) 


To  determine  Q,  two  theorems  must  first  be  proved  for  the  functions  G„,  like  those 
for  Bessel’s  function. 

Theorem  I.  Let  G/  and  Gj  be  two  functions  satisfying 


1  dG 


JT  +  (  i-r*2  +  ^  Q 


(D.3a) 


d^Gj  j  dGj 

^  ^  7*dF*  ^  ^  )G.  =  0 


if  A/  and  Ay  are  the  roots  of  equation  G(l)  =  0,  and  A/  7^  Ay,  then 


G,Gy  r*(  1-r*^  +  AKn  )  dr*  =  0 


(D.3b) 


(D.4) 


Let  A,-  and  Ay  be  the  first  two  eigenvalues.  Multiplying  Eq.  (D.3a)  with  Gjdr*,  Eq. 
(D.3b)  with  Gidr*  and  integrating  from  0  to  1,  and  subtracting  each  other,  we  obtain 

(A,^-A/)  I  G/Gj  r*(  ]-r»2  +  4i:„  )  dr*=  (  G  ^-G-^  ^ 

Jo  ^  'dr*  ^Idr*  ''*-1  (D.5) 


Because  both  G,-  (1)  =  0  and  Gj  (I)  =  0,  A.-  #A,.,  therefore,  Eq.  (D.4)  must  be  true. 
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Theorem  n.  The  value  of  the  integral  in  Eq.  (D.4)  for  i  =  j  can  be  determined  by 


I  ^  Gfi,  r.(  +  4Kn  )  dr*= 


(D.6) 


Let  k  and  Xj  be  the  first  two  eigenvalues  again,  and  the  Eq.  (D.5)  is  true.  Now  let 


A;  -  A;  +  dXi  , 


2A,-  [  GjGi  r*i  l-r*2  +  )  dr*=  _  G— ^ 

J  0  ^  dXi  dr*  ^‘dXidr 

Since  the  special  A,  is  one  of  the  roots  of  equation  G  (1)  =  0,  we  have 

I  Cfi,  r*(  l-r*2  +  4Kn  )  dr»=  4-(  ^  ) 

Jo  2A,- ^  dXidr*  ^^*  =  1 

With  the  help  of  these  two  theorems,  multiplying  the  Eq.  (D.2) 

00 

1  =  2  C„Gn 

n  =  l 

with  Gi  r*(l-  r*2 )  dr*  and  integrating,  we  obtain 


)  r*  =  i 


(D.7) 


r*(  l-r *^  +  4Kn  )  dr* 


=  C,  f  '  Gflj 

J  0 


(D.8) 


.  r*(  l-r*2  +  AKn  )  dr* 


(D.9) 


From  the  differential  equation  (Eq.  (D.3a)),  the  left  side  of  the  integration  equals : 


\2^  dr*  ^  i*  =  \ 
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